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ABSTRACT 


The  objective  of  this  research  is  to  develop  a  mathematical  model  and  control 
algorithm  for  the  maintenance  of  the  environmental  system  within  the  Controlled 
Environment  Research  Chamber  (CERC)  located  at  the  National  Aeronautics  and 
Space  Administration  (NASA)  Ames  Research  Center,  Moffet  Field,  CA.  The 
hypobaric  research  chamber  is  currently  undergoing  renovation  as  part  of  the 
Human  Exploration  Development  Project  (HEDP),  an  effort  on  behalf  of  NASA 
for  advanced  life  support  research.  A  broad  overview  of  the  chamber  is  provided 
which  includes  a  physical  description,  preliminary  system  hardware  and  associated 
performance,  and  potential  experimental  uses.  A  mathematical  model  of  the 
chamber  air  mass  has  been  developed  based  on  key  energy  and  mass  balances. 
Two  methods  of  adaptive  control  have  been  implemented  for  the  coupled  control 
of  temperature,  oxygen  concentration,  pressure  and  humidity  within  the  closed 
environment.  Simulations  testing  algorithm  performance  have  been  conducted, 
including  a  step  and  modified  ramp  response.  The  results  of  the  simulations 
indicate  the  adaptive  methods  performed  well  for  the  model  presented.  Further 
research  is  required  in  refining  the  chamber  model  for  algorithm  optimization  and 
validation  including  the  integration  of  selected  hardware  dynamics. 
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I.  INTRODUCTION 


A.  HUMAN  EXPLORATION  DEMONSTRATION  PROJECT 

The  exploration  of  new  frontiers  has  always  been  a  driving  force  in  human 
evolution.  A  return  to  the  lunar  landscape  and  on  to  the  planet  Mars  appears 
to  be  the  next  logical  step  beyond  the  establishment  of  the  Space  Station 
Freedom.  The  National  Aeronautics  and  Space  Administration  (NASA)  Ames 
Research  Center,  in  support  of  the  Space  Exploration  Initiative  (SEI).  has 
undertaken  an  internal,  multi-division  project  referred  to  as  the  Human 
Exploration  Demonstration  Project  (HEDP).  This  project  has  been  portrayed 
as  an  avenue  to  address  the  advanced  technology  requirements  necessary  to 
establish  an  efficient,  integrated  human  living  environment  on  either  a  lunar  or 
Martian  planetary  surface.  Specifically,  this  high-tech  environment  will  be 
capable  of  providing  all  the  necessary  life  support  functions  in  an  efficient  closed 
loop  system.  Included  within  the  environment  will  be  the  extensive  utilization 
of  state  of  the  art  automation  for  both  physiological  and  psychological 
monitoring,  data  acquisition,  and  overall  habitat  regulating  systems.  HEDP  will 
provide  a  "test  bed"  for  the  research  and  development  of  these  critical 


technologies. 


Four  basic  goals  have  been  identified  for  the  Human  Exploration 
Demonstration  Project:  [Ref.  1.  p.  2]. 

•  Provide  a  simulator  for  the  demonstration  and  evalution  of  selected 
technologies  in  an  integrated  setting. 

•  Create  a  realistic  environment  for  introduction  of  new  technology. 

•  Enhance  the  technology  development  and  evaluation  process  through 
synergistic  cooperation  of  multiple  Ames  divisions. 

•  Identify  promising  technology  concepts  to  programmatic  Centers  for  new 
and  existing  NASA  projects. 

The  spectrum  and  variety  of  potential  experiments  separates  the  HEDP 
project  from  similar  individual  development  and  demonstration  efforts.  Although 
the  central  theme  will  concentrate  on  manned  space  exploration  and  the 
technologies  involved  in  closed  habitats,  the  opportunities  to  develop  systems 
encompassing  the  full  space  mission  spectrum  are  numerous.  The  Human 
Exploration  Demonstration  Project,  in  an  effort  to  meet  all  of  the  aforementioned 
requirements,  will  be  comprised  of  four  major  facilities  and/or  dedicated  systems. 
These  sytems  are  the  Controlled  Environment  Research  Chamber  (CERC).  a 
lunar  landscape  physical  simulation,  robotic  platforms  and  devices,  and  a  human 
powered  centrifuge.  The  majority  of  these  facilities  will  be  contained  within 
existing  buildings  and/or  entities  currently  available  within  the  confines  of  the 
NASA  Ames  Research  Center  Complex,  Moffet  Field,  CA. 


B.  CONTROLLED  ENVIRONMENT  RESEARCH  CHAMBER 

A  critical  component  in  the  overall  HEDP  project  is  the  Controlled 
Environment  Research  Chamber.  Designed  as  a  hypobaric  facilitv'  for  use  in 
high  altitude,  human-based  research  in  partial  support  of  the  NASA’s  drive  to  the 
moon,  it  apparently  went  into  a  state  of  disuse  in  1975-76  [Ref.  1,  p  5].  The 
existing  facility  is  currendy  housed  in  Building  N-239A.  It  consists  of  a  main 
chamber  vessel  and  an  accompanying  airlock  (Figure  1. 1).  The  primary  chamber 


Figure  1.1  Controlled  Environment  Research  Chamber 
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facility  is  a  vertically  oriented  cylinder  with  an  approximate  diameter  of  16  feet. 
The  main  body  extends  two  levels  in  building  N-239A  and  has  an  internal  tloor 
plan  which  consists  of  two  separate  deck  spaces.  The  internal  floorspace 
encompasses  aoout  165  square  feet.  Movement  between  the  levels  is 
accomplished  via  a  wall  mounted  ladder  which  extends  to  the  second  floor.  The 
airlock  was  designed  to  provide  a  buffer  for  entrance  into  the  main  chamber. 
The  airlock  is  a  horizontal  cylinder  that  is  approximately  8  feet  in  diameter  and 
provides  air  tight  entries  from  both  uie  outside  into  the  airlock  and  from  the 
airlock  into  the  main  chamber. 

The  habitat  simulation  will  be  housed  within  the  Controlled  Environment 
Research  Chamber.  The  internal  environment  will  provide  a  test  bed  for  the 
evaluation  of  the  various  subsytems  to  sustain  life  support,  verify  algorithms  and 
methods  for  the  monitoring  of  crew  health,  and  for  the  design  of  elements  for  the 
maintenance  of  personal  hygeine.  A  typical  lunar  and/or  Mars  based  habitat 
must  be  capable  of  providing  certain  minimum  functions.  Those  functions 
include:  [Ref.  1  p.  5]. 

•  Temperature  and  humidity  control 

•  Atmospheric  monitoring  and  compositon  control 

•  Air  revitlization 

•  Water  reclamation 


4 


•  Solid  Waste  Management 

The  HEDP  program  definition  currently  only  addresses  the  first  three  ot  the 
aforementioned  functions.  The  temperature  and  humidity  control  susbsytem  will 
basically  function  by  either  heating  or  cooling  the  constant  ventilation  flow¬ 
through  the  chamber  via  a  counter  flow  condensing  heat  exhanger.  The  humidity 
of  the  air  will  be  controlled  utilizing  the  same  apparams.  The  compostion  ot  the 
air  will  be  monitored  and  the  oxygen-nitrogen  levels  of  the  environment  will  be 
tracked.  Makeup  of  flows  of  oxygen  and  nitrogen  will  be  provided  via  a  supply 
of  cryogenic  gases.  A  system  to  regulate  the  removal  and  control  of  trace 
contaminants,  including  carbon  dioxide,  is  also  planned. 

The  inherent  costs  involved  in  the  transport  of  materials  into  space  demands 
design  and  implementation  of  these  ftmctions  in  the  most  proficient  means 
achieveable.  The  testing  and  evaluation  of  the  current  life  support  technologies 
is  best  suited  to  a  controlled  environment  that  closely  emulates  the  conditions 
under  which  the  system  would  operate.  The  perilous  environment  found  on  both 
the  lunar  surface  and  the  Martian  landscape  dictate  tight,  autonomous  control  of 
the  habitat  atmosphere  system.  The  Controlled  Environent  Research  Chamber 
will  provide  the  avenue  under  which  these  functions  and  control  algorithms  may 
be  demonstrated  and  operated. 
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C.  FUTURE  EXPERIMENTS 


Future  experiments  for  the  Controlled  Environment  Research  Chamber 
consist  of  four  basic  categories.  The  categories  are: 

•  Loop  closure  maximization 

•  System  integration  issues 

•  Impact  of  spacecraft/habitat  design  issues  on  life  support  systems 

•  Bioregenerative  or  hybrid  life  support  systems 

Loop  closure  maximization  will  concentrate  on  the  study  of  new  and  innovative 
life  support  process  techologies  involving  air,  water  and  waste  reclamation  and 
management.  Experiments  involving  systems  analysis  and  integration  will 
revolve  around  the  development  and  evaluation  of  sensors  and  control  hardware, 
the  variation  in  metabolic  loads  on  life  support  overall  system  performance  and 
integration,  and  on  the  long  term  logistic  and  reliability  requirements  of  a  closed 
life  support  system. 

The  study  of  the  variation  in  environmental  parameters  on  the  life  suppon 
system,  to  include  the  effects  of  pressure,  composition,  and  humidity  will 
emphasize  another  category  of  experiments.  In  addition,  the  monitoring  of  trace 
contaminants  and  quality  of  recycled  water  may  be  accomplished  in  this 
controlled  environment.  Finally,  an  examination  of  bioregenerative  systems, 
which  include  the  use  of  plants,  may  be  accomplished.  This  will  provide  an 
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avenue  for  studies  in  plant  growth  optimization  and  determination  of  the  life 
support  system  functions  that  may  be  supported  by  bioregenerative  components. 

D.  PROBLEM  STATEMENT 

The  goal  of  this  thesis  is  to  develop  a  mathematical  model  and  control 
algorithm  for  the  integrated  control  of  the  temperature,  humidity,  pressure,  and 
composition  of  the  habitat  enclosed  within  the  Controlled  Environment  Research 
Chamber.  The  model  will  include  a  basic  review  of  hardware  and  sensors  that 
have  already  been  targeted  for  use  in  the  chamber  overhaul.  A  review  of  the 
expected  oxygen  and  nitrogen  consumption/makeup  requirements  in  addition  to 
an  estimate  of  both  the  sensible  and  latent  heat  loads  will  be  accomplished  and 
a  preliminary  design  to  meet  these  requirements  will  be  proposed.  Concurrently, 
the  aspects  of  oxygen  deprivation  and/or  the  bends  phenomenon  will  be 
researched  to  ensure  the  control  system  is  capable  of  safely  supporting  human 
inhabitants.  Two  different  modem  control  methodologies  will  be  examined  to 
include  an  adaptive  scheme  of  parameter  estimation  and  control  gain  calculation 
and  one  which  includes  the  use  of  an  online  parameter  estimator  coupled  with  a 
linear  quadratic  regulator  control  algorithm. 

The  system  must  be  capable  of  providing  a  suitable  environment  through  a 
varied  set  of  operational  conditions.  In  addition,  the  accuracy  of  control  must 
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be  maximized  to  the  largest  extent  possible  to  fully  utilize  the  chamber  for  a  wide 
tightness  of  control  requirements.  A  summary  is  provided  in  Table  I. 


Table  I  -  CERC  Operating  Regime  and  Control  Requirements 


Parameter 

Nominal 

Operating 

Range 

Absolute 

Min/ Max 

Values 

Nominal  Value 

Control 

Accuracy 

Required 

Pressure 

34.464  - 
101.325  kPa 

34/102  kPa 

Depends  on 

Experimental 

Requirements 

1%  of 
Nominal 

Pressure 

Temperature 

18.33-26.67  "C 

15.55/29.44  “C 

21.1  "C 

Depends  on 

Experimental 

Requirements 

Mass  Fraction 
02 

{Dry  Air) 

0.64903- 

0.2321 

0.649-0.232 

Depends  on 
Operational 
Pressure 

-1-/-  2%  of 
Nominal 

Pressure 

Absolute 

Humidity 

0.0079504  - 
0.0232739 
kg  H.O/kg  dry  air 

25/70% 

Relative 

Humidity 

Nominally 
Controlled  (0 
50%  RH 

Depends  on 

Experimental 

Requirements 
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Critical  to  the  overall  control  scheme  is  maintenance  of  the  partial  pressure 
of  02  when  transitioning  between  experimental  set  points.  Accuracy  of  control 
of  the  temperature  and  humidity  will  be  nominally  set  at  ±  1  '  and  ±5%  of 
absolute  humidity,  respectively. 

The  response  of  the  control  system  will  then  be  analyzed  with  various 
distrubances  and  metabolic  loads.  The  feasiblity  of  an  adaptive  control  algorithm 
will  then  be  reviewed  for  utilization  in  the  future  chamber  design  implementation. 
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11.  BACKGROUND 


The  following  sections  detail  two  major  subject  areas  requiring  discussion 
with  regards  to  the  control  of  the  environmental  conditions  within  the  Controlled 
Environment  Research  Chamber.  Those  areas  include  a  brief  synopsis  on 
humidity  and  some  insight  into  the  impact  of  working  in  a  reduced  pressure 
environment. 

A.  PSYCHROMETRICS 

Psychrometrics  is  the  study  of  the  impact  of  moisture  on  the  properties  of 
air.  The  simplest  manner  to  obtain  an  inference  into  the  moisture  content  of  air 
is  through  an  application  of  Dalton’s  law  of  partial  pressures.  This  law  simply 
states  that  in  any  nonreacting  mixture  of  gases  and  vapors,  each  gaseous  or  vapor 
component  exerts  an  individual  partial  pressure  that  is  equal  to  the  pressure  that 
the  gas  would  exert  if  it  occupied  the  space  alone,  and  that  the  total  pressure  is 
the  sum  of  the  mixture  exerted  by  the  individual  gases  or  vapors  [Ref  2,  p.  21]. 
Hence,  the  total  measured  barometric  pressure  has  a  component  representing  the 
water  vapor  and  is  the  basis  for  the  definition  of  humidity. 
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Ii  is  important  to  recognize  that  the  water  vapor  in  air  is  actually  steam  at 
low  pressure.  The  maximum  amount  of  vapor  that  can  be  mixed  with  any  given 
volume  of  dry  air  depends  only  on  the  temperature  of  the  air.  Since  the  amount 
of  water  vapor  in  the  air  determines  the  partial  pressure  exened  by  the  water 
vapor,  it  is  evident  that  the  air  will  contain  the  maximum  amount  of  water  vapor 
when  the  water  vapor  in  the  air  exerts  the  maximum  possible  pressure.  The 
maximum  pressure  that  can  be  exerted  by  any  vapor  is  known  as  the  saturation 
pressure.  The  saturation  pressure  of  water  vapor  at  any  given  temperature  may 
be  determined  utilizing  the  following  relationship:  [Ref  3,  p.  59] 

logp=28.590l  -8.2log(r+273)+0.00248(7+273) 

where  7 is  the  saturation  temperature  in  X  and  p  is  saturation  pressure  in  bars. 

The  temperature  corresponding  to  saturation  is  commonly  referred  to  as  the 
dew  point  temperature.  The  presence  of  ’dew’  in  the  morning  is  an  indication 
that  the  temperature  dropped  below  the  dew  point  temperature. 

The  water  vapor  content  in  the  air  is  called  humidity.  The  terms  absolute 
and  relative  humidities  are  conventionally  utilized  to  define  the  moisture  state  of 
the  air.  Relative  humidity  (0)  is  the  ratio  of  the  mole  fraction  (or  partial 
pressure)  of  water  in  moist  air  to  the  mole  fraction  (or  partial  pressure)  of  water 
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in  saturated  air  at  the  same  temperature  and  pressure.  Relative  humidity  (RH) 
can  be  expressed  functionally  as 


(2.2) 


where  p^.  is  the  measured  partial  pressure  of  the  water  in  the  air 

is  the  saturation  partial  pressure  of  the  water  in  the  air 
Absolute  humidity,  on  the  other  hand,  is  defined  at  the  acmal  density  of 
water  vapor  in  air.  For  many  applications,  it  simplifies  the  analysis  to  define 
absolute  humidity  using  a  term  known  as  the  humidity  ratio,  w.  This  is  simply 
the  ratio  of  the  mass  of  water  to  the  mass  of  dry  air.  The  humidity  ratio  at 
saturation  may  by  calculated  via  the  following  relationship  [Ref.  3,  p.20]. 


K  P.rPs 


(2.3) 


where  is  the  molecular  weight  of  H2O  (kg/kg-mole) 

is  the  molecular  weight  of  dry  air  (kg/kg-mole) 
p,  is  the  partial  pressure  of  H2O  in  air  (kPa) 

Pa  is  the  total  atmospheric  pressure  of  air  (kPa) 

Henceforth,  any  reference  to  absolute  humidity  will  be  actually  a  reference  to  the 
humidity  ratio. 
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The  methods  of  humidity  control  vary.  The  most  common  method  of 
dehumidification  is  the  utilization  of  an  open  coil  condensing  heat  exchanger. 
Moist  air  is  made  to  come  in  contact  with  a  cooling  coil,  which  is  held  at  a 
temperature  below  the  dew  point  of  the  air.  The  moist  air  is  cooled  and  becomes 
saturated.  The  amount  of  moisture  removed  via  condensation  is  a  function  of  the 
mean  coil  temperature  and  the  air  dew  point  temperature.  The  difficulty 
encountered  in  dehumidification  lies  in  the  coupling  of  humidity  and  temperature. 
This  is  because  a  very  small  change  in  temperature  can  cause  a  very  large  change 
in  humidity.  The  easiest  manner  in  which  to  conceptualize  the  dehumidification 
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process  is  to  follow  it  on  the  psychrometric  chan  as  depicted  in  Figure  2.1.  .411 
processes  on  the  chart  are  at  constant  pressure.  Starting  at  a  point  on  the 
saturation  (100%  RH)  curve,  relative  humidity  varies  from  left  to  right.  An 
increase  in  the  sensible  heat  load  is  equivalent  to  a  horizontal  line  in  the  same 
direction.  It  can  be  seen  that  starting  from  a  point  on  the  saturation  curve, 
raising  the  temperature  results  in  lower  humidity.  Moreover,  it  can  be  seen  that 
the  relationship  between  dry  bulb  temperature  and  relative  humidity  is 
exponential.  This  accounts  for  the  fact  that  any  change  in  temperature  results  in 
a  much  larger  change  in  humidity. 

B.  HYPOBARIC  CHAMBER  DESIGN 

Large  decompression  chambers  and  simulated  habitats  present  a  number  of 
potential  hazards  to  the  personnel  enclosed  within  them.  From  an  atmospheric 
control  point  of  view,  one  major  consideration  is  the  decompression  effects 
derived  from  the  improper  transition  to  lower  pressures  and  the  impact  of  oxygen 
deprivation  (hypoxia)  or  excess  (hyperoxia).  An  artificial  atmosphere  of  suitable 
composition  and  pressure  is  critical  to  the  health  and  performance  of  crew 
members.  This  atmosphere  supplies  the  oxygen  their  blood  must  absorb  and  the 
pressure  their  body  fluids  require. 
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Humans  are  accustomed  to  living  in  an  environment  that  contains  21% 
oxygen  by  volume  at  14.7  psia  ambient  pressure.  According  to  Dr.  R.  Vann. 
Duke  University,  the  optimum  oxygen  content  is  obtained  when  the  partial 
pressure  of  oxygen  is  maintained  at  a  nominal  value  of  160  mm  Hg  for  all 
operating  pressures  below  standard  atmospheric  pressure  [Ref.  4],  Large 
variances  from  this  value  can  cause  physiological  complications  such  as  hypoxia 
or  hyperoxia.  Too  little  oxygen  (hypoxia)  induces  sleeplessness,  headaches,  the 
inability  to  perform  simple  tasks,  and  eventually  unconsciousness.  Acute 
impairment  of  brain  function  occurs  within  13  seconds  whenever  the  aveolar 
oxygen  tension  drops  below  33  mm  Hg.  Too  much  oxygen  (hyperoxia)  can  lead 
to  respiratory  problems  including  inflammation  of  the  lungs,  heart  disturbances, 
blindness,  and  loss  of  consciousness  [Ref.  5,  p.  5-1].  Figure  2.2  shows  human 
performance  limits  versus  total  pressure  and  oxygen  concentration. 

In  addition,  the  maintenance  of  the  pressure  and  the  proper  timing  of 
pressure  transitions  is  crucial.  There  must  exist  a  total  pressure  that  prevents  the 
vaporization  of  body  fluids.  This  is  referred  to  as  an  ebulism  and  occurs  around 
one  psia  at  a  temperature  of  37°C  [Ref.  5,  p,  5-4].  The  transition  from  an 
atmospheric  state  to  a  hypobaric  state  in  an  improper  or  excessivley  rapid  manner 
can  cause  decompression  sickness.  Decompression  sickness  is  the  formation  of 
small  gas  bubbles  in  the  body  tissues  or  vascular  system  and  is  commonly 
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Figure  2.2  -  Oxygen/Pressure  Envelopes  for  Optimum  Human  Performance 


[Ref.  5,  p.  5-1] 


referred  to  as  the  bends.  The  seriousness  of  a  case  of  the  bends  is  related  to  an 


anatomic  lodgement  of  one  of  these  bubbles  which  may  press  on  nerves  or 
obstruct  the  blood  supply.  NASA  currently  utilizes  a  decompression  scenario  for 
Extra  Vehicular  Activities  (EVA)  in  space  suits  during  shuttle  missions  which 
requires  a  staged  decompression.  The  pressure  is  reduced  to  around  10.2  psia 
where  a  period  of  pure  oxygen  pre-breath  is  conducted.  The  transition  to  the 
space  suit  pressure  is  then  accomplished  [Ref  6,  p.  5  ]. 

For  the  likely  scenarios  involving  the  Controlled  Environment 
Research  Chamber,  this  was  taken  as  the  standard  decompression  scheme.  A 
candidate  control  system  must  thus  be  capable  of  accurately  meeting  this 
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requirement  while  maintaining  the  oxygen  partial  pressure  at  a  nominal  value  of 
160  mm  Hg. 

From  a  controls  point  of  view,  it  is  imperative  to  accurately  and 
properly  transition  from  one  pressure  state  to  another.  In  addition,  the 
concentration  of  oxygen  must  be  maintained  precisely  to  avoid  either  oxygen 
deprivation  (hypoxia)  or  oxygen  excess  (hyperoxia).  The  wide  range  of 
operating  pressures  and  varied  experimental  requirements  for  the  Controlled 
Environment  Research  Chamber  further  complicates  the  control  algorithm 
inception. 
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III.  MATHEMATICAL  MODELING 


A.  SYSTEMS  OVERVIEW 

The  atmospheric  control  system  is  comprised  of  two  distinct  subsystems. 
One  of  the  subsystems  will  provide  the  monitoring  and  control  of  the  ambient 
pressure  and  concentration  of  oxygen.  A  cryogenic  gas  supply  system  will 
provide  makeup  gas  with  a  vacuum  line  assisting  in  pressure  control.  This  is 
known  as  the  Atmospheric  Monitoring  and  Control  Subsystem.  The  temperature 
and  humidity  will  be  maintained  via  the  Temperature  and  Humidity  Control 
Subsystem.  This  system  consists  of  temperature  and  humidity  sensors  coupled 
with  a  condensing  heat  exchanger  capable  of  simultaneously  removing  excess 
moisture  and  maintaining  chamber  temperature. 

1.  Atmospheric  Monitoring  and  Control  Subsystem 

The  Atmospheric  Monitoring  and  Control  Subsystem  is  specifically 
tasked  with  the  maintenance  of  the  proper  mix  of  atmoshperic  gaseous 
constituents  within  the  limits  set  for  experimental  requirements.  In  addition,  the 
system  is  responsible  for  steady  state  pressure  and  transitional  pressure  control. 
The  system  is  to  provide  autonomous  operations  utilizing  a  computer  based 
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system  with  the  capability  to  manually  control  the  system  in  the  event  of  an 
emergency  condition. 

The  monitoring  of  the  gaseous  components  will  be  accomplished  via 
an  Ohmeda  Rascal  II  Anesthetic  Gas  Monitor.  The  Rascal  II  uses  laser  Raman 
scattering  characteristics  of  gases  to  identify  the  gaseous  constiments  and  relative 
concentrations.  This  device  is  capable  of  analyzing  up  to  eight  atmospheric  gases 
with  a  minimum  detectability  of  0.25%  by  volume  at  a  time  response  of  less  than 
or  equal  to  500  msec  [Ref.  7,  p.4].  Multiple  sensors  may  be  employed 
throughout  the  chamber  and  read  through  a  multiplexer  to  obtain  a  more  accurate 
overall  view  of  the  gaseous  distribution.  Concentrated  pockets  of  a  particular 
constituent  may  be  detected  utilizing  this  methodology  and  a  potentially 
hazardous  situation  avoided.  A  pressure  sensor  has  not  been  chosen  for  the 
chamber  renovation. 

Make  up  gases  will  be  provided  via  a  Linde  cryogenic  gas  supply 
system.  The  gases  will  be  mixed  in  a  gas  blending  unit  with  the  resulting  mix 
injected  into  the  chamber  via  the  supply  ducting  in  the  ventilation  system. 
Control  of  the  various  flows  will  be  accomplished  via  a  combination  of  a  Linde 
flow  meter  console  and  Linde  flow  control  modules.  Each  operator/flow 
metering  console  is  capable  of  controlling  the  flow  rates  in  four  separate  control 
modules.  The  operator  console  is  configured  with  a  RS232  connection  for 
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external  operation  from  a  microprocessor  or  hands-off  programmed  operation  by 
a  central  computing  station.  Manually  operated  valves  will  be  provided  for 
emergency  operations.  In  addition,  a  completely  separate  oxygen  supply  system 
consisting  of  an  independent  piping  system  and  gas  masks  available  within  the 
confines  of  the  chamber  will  be  provided  as  an  emergency  breathing  alternative. 


2.  Temperature  and  Humidity  Control  Subsystem 

The  Temperature  and  Humidity  Control  Subsystem  is  designed  to 
enhance  human  comfort  by  controlling  the  gas  temperature  and  absolute 
humidity.  It  is  also  configured  for  autonomous  operations  with  the  capability  to 
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be  maunually  overriden  in  the  event  of  an  emergency.  A  temperature  sensor  has 
not  been  chosen  for  the  chamber  overhaul  process.  For  modeling  purposes,  a 
thermocouple  has  been  selected  for  the  chamber  renovation.  A  thermocouple  has 
an  extremely  fast  response  time  ( <  1  sec)  and  requires  no  external  power 
source.  This  fast  response  time  will  alleviate  evenmal  modeling  of  a  sensor  time 
delay.  The  signal  does,  however,  require  some  amplification  and  the  output 
would  require  some  additional  conditioning  due  to  nonlinearities  in  the  output. 
For  the  tracking  of  humidity,  a  modem  humidity  sensor  will  be  employed. 
Modem  humidity  sensors  employ  thin  film  technology  and  measure  the  ambient 
air  dew  point  (from  which  the  partial  pressure  may  be  derived)  as  a  change  in  the 
capacitance  of  either  aluminum  hydroxide  or  silicon. 

The  ventilation  system  consists  of  an  eight  inch  exhaust  line  and  a 
twelve  inch  supply  line.  The  internal  ducting  will  be  sized  and  configured  to 
ensure  that  the  air  flow  does  not  interfere  with  normal  experimental  procedures. 
The  air  mover  is  a  sealed  centrifugal  fan  which  will  have  a  nominal  capacity  of 
3400  cubic  feet  per  minute.  The  air  flow  will  be  channeled  through  air  tight 
ducting  into  a  sealed  condensing  heat  exchanger.  A  damper  mechanism  will  be 
utilized  to  control  or  divert  the  required  flow  over  the  condensing  heat 
exchanger.  The  diverted  flow  will  then  be  remixed  with  the  heat  exchanger  flow 
prior  to  injection  back  into  the  chamber.  The  temperature  will  be  primarily 


21 


controlled  via  the  magnitude  of  air  flow  through  the  condensing  heat  exchanger. 
The  humidity  will  be  primarily  controlled  via  coolant  entrance  temperature. 
Figure  3.2  provides  a  schematic  of  the  Temperature  and  Humidity  Control 
subsystem. 


Figure  3.2  Temperature  and  Humidity  Control  Subsystem 


A  specific  condensing  heat  exchanger  has  not  been  chosen  for  the 
chamber  renovation  project.  However,  parametric  studies  were  conducted 
utilizing  a  conventional  cross  flow  heat  exchanger  with  a  nominal  heat  transfer 
coefficient  (U)  of  10  btu/hr'ft^-°F.  Sensible  and  latent  heat  loads  utilized  in  the 
studies  were  derived  from  information  contained  in  reference  6  and  will  be 
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detailed  in  the  following  section.  A  25%  solution  of  calcium  chloride  brine  was 
selected  as  the  heat  exchanger  coolant  because  of  its  low  freezing  point  (-  1 1  'C) 
and  non-toxic  nature  [ref.  9,  p.  5].  The  performance  of  the  heat  exchanger  was 
defined  in  terms  of  nominal  coolant  rise  temperature  and  median  approach 
temperature.  The  nominal  coolant  rise  temperature  was  modeled  at  5.55  'C. 
Approach  temperature,  the  difference  between  the  temperamre  of  the  air  that 
exits  the  heat  exchanger  and  the  coolant  entrance  temperamre,  was  set  at  8.33 
°C.  The  parametric  smdies  were  conducted  at  various  pressures  between  34.4 
and  101.325  kPa  and  the  results  are  provided  in  Appendix  B  [Ref.  10].  The 
required  heat  exchanger  was  nominally  sized  at  5.6  m^  of  frontal  area  with  a 
coolant  flow  rate  of  178  kg/hr.  These  parameters  were  utilized  in  the 
mathematical  modeling  of  the  overall  atmospheric  control  system. 

B.  DYNAMIC  EQUATIONS  DEVELOPMENT 

The  state  of  the  air  mass  in  the  chamber  at  any  given  instant,  given  the 
assumptions  provided  in  the  following  section,  can  be  fully  specified  with  four 
state  variables.  Those  variables  are  pressure,  mass  fraction  of  oxygen  in  dry  air. 
temperamre  and  absolute  humidity.  The  following  subsections  will  detail  the 
assumptions,  nomenclamre  and  equations  utilized  in  the  development  of  the 
chamber  model  representing  the  dynamics  of  these  four  variables. 
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1.  Assumptions 


The  following  is  a  listing  of  the  major  assumptions  made  in  the 
development  of  the  equations  defining  the  air  mass  dynamics  within  the  Closed 
Environment  Research  Chamber. 


1.  The  gaseous  consituents  in  the  chamber  are  assumed  to  behave  as  ideal 
gases. 

2.  The  chamber  air  is  composed  of  only  oxygen,  nitrogen  and  water  vapor. 
All  contaminants  and  minor  gases  are  excluded  from  consideration. 

3.  Potential  energy  of  the  chamber  gas  mass  is  assumed  to  be  negligible. 

4.  The  chamber  is  modeled  as  a  lumped  system  to  simplify  analysis  by 
avoiding  partial  derivatives. 

5.  The  performance  of  the  heat  exchanger  is  fixed  and  is  not  varying  over 
time. 

6.  The  chamber  is  housed  within  a  climate  controlled  building  and  hence  the 
heat  transfer  through  the  walls  into  or  out  of  the  chamber  is  considered 
negligible. 

7.  The  interior  mean  radiant  and  structure  temperatures  are  assumed  to  be 
equal. 

8.  There  are  no  objects  in  the  chamber  that  absorb  heat  or  moisture. 

9.  The  velocity  of  the  airflow  through  the  chamber  is  assumed  to  be  in  a 
uniform  direction  and  does  not  appreciably  affect  pressure  (static 
pressure = stagnation  pressure) . 

10.  Ventilation  flow  through  the  chamber  is  assumed  constant  at  3400  cubic 
feet  per  minumte  (cfin). 
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11.  It  is  assumed  that  all  makeup  gases  enter  the  chamber  at  ambient 
temperature  and  pressure  and  are  perfectly  mixed. 

12.  Particle  concentrations,  pressure,  humidity  and  temperature  are  measured 
at  the  perfectly  mixed  conditions. 

13.  The  heat  capacities  and  thermal  properties  of  the  chamber  constituents  are 
assumed  to  be  not  a  function  of  pressure. 

14.  The  production  of  CO,,  derived  through  the  metabolism  of  O,,  is  assumed 
to  be  removed  by  the  COi  removal  subsytem  as  it  is  produced  and  hence  is 
neglected. 

15.  It  is  assumed  that  there  is  no  delay  in  the  actuation  of  the  controls  and  that 
the  inputs  have  an  instantaneous  effect  on  the  chamber  states. 

16.  The  maximum  makeup  flow  rates  of  oxygen  and  nitrogen  and  vacuum  are 
assumed  to  be  10  kg/hr. 

17.  The  temperature  range  of  coolant  is  assumed  to  be  between  0  and  15  °C. 

18.  The  chamber  is  assumed  to  have  four  human  inhabitants  with  the  following 
oxygen  consumption  and  heat  production  loads.  Three  of  the  humans  are 
assumed  to  be  engaged  in  normal  activity  while  one  is  involved  in 
exercise/strenuous  activity  [Ref.  8,  p.  10  ]. 

Table  II  SUMMARY  OF  CHAMBER  LOADS 


Source 

Nominal  O2 
consumption  (kg/hr) 

Sensible  heat 
production  (kJ/hr) 

Latent  heat  (H.O) 
production  (kg/hr) 

Three  persons, 
normal  activity 

0.035 

1326 

0.285 

One  person, 
strenuous  activity 

0.349 

664 

0.455 

Miscellaneous 
(Machinery, 
lighting,  etc) 

0,0 

7720 

0.0 

Total 

0.3848 

9710 

0.74 
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2.  Nomenclature 


The  following  is  a  listing  of  and  definitions  for  the  variables  utilized 
in  the  system  modeling. 

a,,a2,a3,b,,b2,b3=  Coefficients  of  heat  capacities  for  oxygen,  nitrogen  and 
water. 

£5=  Total  chamber  air  mass  energy  (kJ) 

hvap=  Heat  of  vaporization  of  moisture  generated  by  humans  (kJ/kg) 
hvapl=  Heat  of  condensation  of  moisture  in  heat  exchanger  (kJ/kg) 

Hv=  Heat  of  vaporization  at  0°C  (kJ/kg) 

H20=  Total  human  latent  load  (moisture)  (kg/hr) 
nicHx=  Mass  flow  rate  of  air  through  heat  exchanger  (kg/hr) 
iflLEAK=  Leakage  into  the  chamber  (kg/hr) 
mN2=  Mass  flow  rate  of  nitrogen  into  chamber  (kg/hr) 
mo2=  Mass  flow  rate  of  oxygen  into  chamber  (kg/hr) 

•Mout=  Mass  flow  rate  of  vacuum  flow  out  of  chamber  (kg/hr) 

M,=  Total  chamber  dry  air  mass  (kg) 

Mc=  Total  chamber  air  mass  (kg) 

Mo2=  Total  mass  of  oxygen  (kg) 

M^=  Total  mass  of  rnoismre  (kg) 

02=  Total  human  consumption  of  oxygen  (kg/hr) 
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P,  =  Chamber  Pressure  (kPa) 

p=^=  Partial  pressure  of  water  at  saturation  (kPa) 

q^=  Miscellaneous  sensible  loads  (kJ/hr) 

qh=  Total  human  sensible  heat  load  (kJ/hr) 

qcHx—  Total  heat  removed  in  heat  exchanger  (kJ/hr) 

R=  Universal  gas  constant  8.314  kJ/kg-moIe-  K 
R’=  Specific  gas  constant  (kj/kg  -  K) 

T(.=  Chamber  temperature  (“C) 

Tj=  Temperature  of  coolant  entering  heat  exchanger  {°C) 
Texii=  Temperature  of  air  exiting  heat  exchanger  (°C) 
u= Control  variables 

Chamber  volume  (m^) 

=  Chamber  absolute  humidity  (kg  HjO/kg  dry  air) 
x=  State  variables 
Xj  =  Mass  fraction  of  Oj  in  dry  air 
1-Xj“  Mass  fraction  of  N2  in  dry  air 
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3.  State  Variable  Definition 
a.  Temperature 

Determination  of  the  chamber  bulk  temperature  involves  an 
application  of  the  generalized  conservation  of  energy  equation  (open  system 
formulation)  in  which  the  kinetic  and  potential  energy  terms  are  neglected. 

Lrhh.=Lnih+—  (3.1) 

'  '  '  W/ 

The  terms  L/Wj/ijand  represent  system  heat  inputs  and  outputs,  respectively, 
and  E  represents  the  total  air  mass  energy.  The  total  energy  of  the  system  is 
defined  as  the  mass  of  air  multiplied  by  its  specific  enthalpy.  The  mass  of  air 
may  be  derived  utilizing  the  perfect  gas  law. 

Py^nRiT^^m)  (3.2) 

Because  the  humidity  and  mass  fraction  of  oxygen  are  all  defined  in  terms  of  dry 
air,  the  determination  of  the  system  energy  will  require  referencing  the  system 
air  mass  to  dry  air.  Rearranging  the  ideal  gas  law  and  inserting  the  specific  gas 
constant  R '  for  the  unversal  gas  constant,  the  chamber  dry  air  mass  is  defined  as 

M  = - _ !—  (3.3) 

'  R'{Tym)  1+vv^ 
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where 


R  (3.4) 

"(x.(32)+(1-jc^)(28)-^w^(18)) 

The  heat  capacity  of  the  overall  air  mass  can  be  expressed  utilizing 
the  heat  capacities  of  the  individual  gaseous  components  and  multiplying  by  their 
respective  mass  fractions. 


Cpair=x^(a^+bJ^)H^-x^)(a^+bJc)+w^{aj*bJ^)  (3.5) 

All  of  the  above  heat  capacities  are  based  on  a  specific  temperature 
range  with  a  minimum  of  O^C  (Ref.  2,  p,4%].  The  enthalpy  of  the  chamber  air 
is  defined  as  the  sum  of  the  heat  capacities  of  the  gaseous  components  multiplied 
by  the  chamber  temperature  and  the  heat  of  vaporization  of  water  based.  The 
heat  of  vaporization  for  this  application  will  be  referred  to  as 


(3.6) 


Combining  equations  3.3  and  3.6,  an  expression  for  the  total  energy  of  the 
system  may  be  obtained. 
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£,= - i_f - L.[x  {a,T^bjhHi-x  ){aj  ^b.T^)^w  {H  ^aj  +bj;)] 

'  R'(r+273)  l+w  "  '  "  '  "  "  -  ^  ^  ' 


(3.7) 


where  H,=  2213.35  kJ/kg  [Ref.  1,  p.  395). 

Equation  3.7  must  be  differentiated  with  respect  to  each  of  the 
system  variables  in  order  to  obtain  an  expression  for  dEJdt.  Utilizing  the  chain 
rule,  equation  3.7  was  differentiated  with  respect  to  P„  T^,  and  The  result 
was  then  combined  with  the  heat  inputs  and  outputs  into  (3.1).  For  this  analysis, 
the  heat  inputs  are  defined  as 


i:h/n.=q^*q„*H^O  *  hvap 


and  the  heat  outputs  are  defined  as 


(3.8) 


(3.9) 


The  term  is  defined  as  the  total  heat  (latent  and  sensible)  removed  from  the 
chamber  air  mass  via  the  condensing  heat  exchanger.  As  presented  previously, 
the  nominal  exit  temperature  of  the  air  leaving  the  heat  exchanger  has  been 
modeled  based  on  the  approach  temperature,  which  was  assumed  to  be  a  measure 
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of  performance  for  the  condensing  heat  exchanger.  Assuming  the  air  exit 
temperature  is  a  function  of  this  approach  temperature,  the  partial  pressure  of  the 
moisture  in  the  exit  stream  can  be  determined  via  equation  (2.1)  (assuming  the 
air  is  samrated).  The  heat  of  vaporization  (condensation)  for  the  exit  stream  can 
then  be  calculated  via  an  application  of  the  Clausius-Clapeyron  equation 

hvap={n.  1 1  -2.3Iog/’  *)/?'(7  +273)  (3-10) 

where  P*  is  in  atmospheres  and  hvap  is  in  kJ/kg  [Ref  1,  p.  15]. 

From  equation  (2.2),  the  exit  humidity  can  be  calculated.  Incorporating  this  exit 
stream  humidity  into  the  following  relationship  provides  the  latent  heat  removed 
across  the  condensing  heat  exchanger. 

q=J2^(w-w)hvap  (3.11) 

where  is  the  exit  humidity.  The  sensible  heat  removed  can  be  expressed  as 
the  difference  of  the  heat  capacities  of  the  entrance  and  exit  streams  multiplied 
by  their  respective  temperatures. 


m 


q  = — —[cp  T  -cp  T  ] 

1  +W  ^  ^ air~ent  c  r  air~cxn  exit* 


(3.12) 
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where  cp^^-em  and  represent  the  entrance  and  exit  heat  capacities,  as 

defined  previously.  Incorporating  these  definitions  into  equation  (3.1),  an 
expression  for  the  dynamics  of  temperamre  with  respect  to  the  other  system 
variables  may  be  obtained. 


m 


i  ‘•‘VV 


R'(.l*w^KT^*n3) 


P  V 

C  C 


I 


1  +w 


UMt  T^*bJ^H\-x^)(a.T.  *bjl)  ■ 


-(a,T^*b,n-a,Tra,Ti]^- 


dw 

C 

~W 


dt 


dt 


IQ 


(3.13) 


where  Q  is  defined  as 


i 
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Q- 


aj  bjl 

a,- — Lj_+2Z?,r-.  '  ‘ 

1  T'  .  1 


r +273 

c 


7+273 


^(1-^ ) 


fl.7  _  bj: 

a.  -  . ^Ib^T,- 


T.+m 


7^+273 


+  H' 


H 


7^+273 


+a  ' — Lj_- — f - 

^  7  +273  7  +273 


(3.14) 


b.  Pressure 

The  pressure  was  modeled  by  invoking  an  overall  mass  balance  on 
the  chamber  atmosphere.  All  mass  consumption  and  production  terms,  including 
both  the  consumption  and  supply  of  oxygen,  supply  of  nitrogen,  air  leakage  into 
the  chamber,  removal  through  the  vacuum  line,  and  the  moisture  generation  and 
removal  terms  were  included.  The  generalized  mass  balance  is  defined  as 


dM^ 

~dr 


m 


(3.15) 


where  and  represent  the  aforementioned  mass  inputs  and  outputs, 
respectively.  Equation  (3.14)  can  be  expanded  to  include  the  respective 
chamber  mass  loads. 
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(3.16) 


dM 

dt 


I  +H' 


(H’ 


wJ-0. 


An  application  of  equation  (3.2).  the  ideal  gas  law,  provides  a  relationship  for 
the  overall  chamber  air  mass. 


P  V 

c  c 

R'(T.^273) 


(3.17) 


Differentiating  equation  (3.16),  utilizing  the  chain  rule,  with  respect  to  pressure 
and  temperature,  and  solving  for  the  differencial  of  pressure  with  respect  to  time 
provides  the  equation  of  motion  for  pressure. 


dP^_R\T^^213) 

_ 


1  ‘ 


dr 


r^+273  dt 


(3.18) 


c.  Mass  Fraction  of  Oxygen 

Determination  of  the  dynamics  of  the  chamber  oxygen  mass 
fraction  involves  another  application  of  the  generalized  conservation  of  mass 
equation. 
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(3.19) 


The  inputs  and  outputs,  m,  and  w„^ ,  refer  to  the  oxygen  generation  and 
consumption  terms.  The  total  oxygen  in  the  air  may  be  expressed  as 


A/_,=x  M  (3.20) 

02  c  a 

where  is  the  total  mass  of  oxygen  and  is  the  mass  of  dry  air  as  defined 
in  equation  (3.3).  Expanding  equation  (3.19)  and  inserting  equation  (3.20),  the 
derivative  of  pressure  becomes 


(3.21) 

Utilizing  the  chain  rule  and  differentiating  with  respect  to  T^,  vv^  and  x^,  and 
solving  for  dx^dt. 


(l+w^)^'(r^+273) 

y~p  ^^02^^P^LEAK  ^c^OUT 

c  c 

+___£ _ 1  + L_ ! 

(r^+273)  dt  1+w^  dt 


x^dP^ 
P,  dt 

c 


(3.22) 
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d.  Humidity 


The  dynamics  of  the  moisture  can  be  determined  via  yet  another 
application  of  a  conservation  of  mass,  this  time  on  the  moisture  component. 


dM 

w 

~dr 


(3.23) 


where  is  the  total  mass  of  the  moisture  in  the  air.  can  be  rewritten  in 
terms  of  the  humidity  ratio  and  mass  of  dry  air  as 

M  =wM  (3.24) 

w  c  a 

The  inputs  into  the  system  are  the  generation  of  moisture  by  humans,  and 
the  output  is  the  moisture  removed  by  the  condensing  heat  exchanger.  These 
are  defined  as 


Incorporating  equation  (3,3),  which  defines  the  mass  of  dry  air,  the  differential 
of  with  respect  to  time  can  be  expanded 
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(3.27) 


dM 

M 

~ir 


wPV 

c  c  c 


(I  (7^273) 


dt 


Taking  the  derivative  with  respect  to  P^,  T^,  and  w^,  and  incorporating  the 
expressions  for  the  moisture  inputs  and  outputs,  the  dynamics  of  the  humidity  can 
be  expressed  as 


1 

II 

/?'(r>273)(Uw; 

H,0 

1  dP  1  dT. 

Pv^w^ 

1  +w 

c 

P^  dt  7^+273  dt 

dt 

r  1  _  1  1 

w^+1 

(3.28) 


4.  State  Space  Representation 

Assimilation  of  the  equations  representing  the  dynamics  of  pressure, 
mass  fraction  of  oxygen,  temperature  and  humidity  into  a  simple  set  of  equations 
was  not  feasible.  The  system  is  highly  nonlinear,  coupled  and  cannot  be  easily 
transposed  into  the  standard  state  space  form 
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x=Ax-^Bu^Cw 


(3.29) 


where  x  is  the  vector  of  state  variables,  u  is  the  vector  of  control  inputs,  and  w 
is  the  disturbance  vector,  which  in  this  cases  is  the  system  loads.  The  state 
vector  will  be  henceforth  defined  as  follows 

(3.30) 


and  the  control  inputs  as 

and  the  disturbance  inputs,  which  are  defined  as  H2O,  O2,,  qt  and  q,. 

For  simulation  purposes,  the  system  was  separated  into  the  following  scheme 

x=F  x+Gix,u)+H{x,w)  (3.32) 

The  matrices  F,  G  and  H  are  presented  in  the  computer  code 
available  in  Appendix  A.  The  matrix  system  was  manipulated  and  solved  for 
dx/dt  in  the  following  manner. 


^OVT 

m, 


'02 
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N2 


m 


CHX 


(3.31) 
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X -Fx =G{x,u)  +H(x,w) 

(3.33) 

x{I-F)=G(x,u)+H(x,w) 

(3.34) 

x=U~F)'HG(x.u)+H(x,w)) 

(3.35) 

The  system  of  equations  were  run  to  observe  the  dynamics  of  the 
modeled  chamber.  The  results  of  the  simulation  were  then  analyzed  to  determine 
the  maximum  integration  period  required  to  ensure  that  the  system  dynamics 
would  be  accurately  represented.  The  integration  period  that  accomplished  this 
goal  was  parametrically  determined  to  be  15  seconds  and  a  control  sampling  time 
of  one  and  one-half  minutes. 
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rv.  CONTROL  LAW  DEVELOPMENT 


A.  CLASSICAL  FEEDBACK  CONTROL 

The  system  of  equations  developed  to  model  the  chamber  air  mass  are  multi- 
variable  and  non-linear  in  nature.  Classical  control  theory  requires  a  multi- 
variable  system  to  be  decomposed  into  single-input/single-output  (SISO)  systems 
with  separate  transfer  functions  linking  each  combination  of  input  and  output 
signals.  In  other  words,  a  system  with  three  inputs  and  three  outputs  has  nine 
transfer  functions  which  must  be  arranged  in  matrix  form.  The  poles  of  the  multi- 
variable  system  are  the  poles  of  all  the  individual  transfer  functions.  The  poles  of 
the  overall  system  is  the  collection  of  poles  of  the  individual  subsystems.  The 
zeros  of  the  multi-variable  system,  however,  do  not  correspond  to  the  zeros  of  the 
individual  transfer  functions.  The  use  of  SISO  techniques  to  design  each 
individual  subsystem  with  the  expectation  that  the  overall  results  will  achieve  the 
desired  level  of  performance  is  unlikely.  Hence,  the  use  of  modern  control 
methods  for  this  project  has  been  pursued. 

B.  ADAPTIVE  CONTROL 

Linear  control  techniques  are  well  defined  and  mature.  Implementation  of 
strict  linear  control  practices  on  non-linear  systems  does  not  always  give 
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satisfactory  results.  The  response  characteristics  of  non-linear  processes  utilizing 
linear  control  systems  may  vary  significantly  because  they  are  normally  linearized 
about  one  operating  condition.  The  need  to  find  a  more  sophisticated  controller 
which  could  automatically  adapt  itself  to  the  changing  characteristics  of  the 
controlled  process  has  lead  to  the  development  of  adaptive  control  routines.  For 
the  application  in  this  project,  the  term  "adaptive  control"  represents  one  particular 
approach  to  designing  a  controller  for  the  non-linear  stochastic  system.  This 
approach  consists  of  assuming  that  the  controlled  process  satisfies  linear  equations 
about  a  particular  operating  condition  only  the  values  of  the  dynamic  coefficients 
defining  the  equations  are  unknown.  It  is  based  on  developing  a  routine  which 
estimates  these  coefficients  and  utilizing  the  resulting  estimates  in  the  desired 
control  law. 

Although  the  system  is  nonlinear,  it  can  be  modeled  using  a  linear  state  space 
representation  about  a  nominal  operating  point.  A  parameter  estimation  scheme 
can  then  be  used  to  identify  the  system  parameters,  assuming  that  the  system  can 
be  considered  linear  around  any  given  operating  point.  In  a  linear  system,  a 
controlled  plant  is  described  by  the  state-vector  differential  equation 

(4.1) 

x=[A]x+{B]u 

y=[C\x^[D]u 
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where 


X 

u 

V 

[A],  [B],  [C],  [DJ 


is  the  state  variable  vector 
is  the  state  input  vector 
is  the  state  output  vector 
are  the  coefficient  matrices 


The  matrices  A,  B,  C  and  D  represent  a  model  of  the  system.  Key  to 
obtaining  an  effective  controller  is  developing  or  estimating  the  values  for  these 
matrices,  which  can  be  unknown  or  time  varying,  depending  upon  the  particular 
operating  condition.  In  an  adaptive  scheme,  the  components  of  the  unknown 
matrices  are  estimated  on-line  by  a  suitable  estimator.  Based  on  the  estimates  of 
these  parameters,  a  control  scheme  such  as  a  linear  quadratic  regulator  may  be 
implemented  for  system  control.  This  class  of  adaptive  control  is  known  as 
indirect  adaptive  control.  Figure  4.1  provides  a  schematic  of  this  method. 

The  direct  mode  of  adaptive  control  is  one  which  attempts  to 
parameterize  the  unknown  system  directly  in  terms  of  the  necessary  control  inputs 
in  order  to  implement  the  desired  control  algorithm.  For  this  mode,  an  estimator 
not  only  provides  system  identification  but  also  predicts  the  control  parameters  and 
utilizes  them  directly  in  the  controller.  The  basic  schematic  is  provided  in  Figure 


4.2. 
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Figure  4.1.  Indirect  Adaptive  Control 


Figure  4.2  Direct  Adaptive  Control 
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C.  INDIRECT  ADAPTIVE  CONTROL 


1.  System  Identification 

Parameter  estimation/system  identification  is  the  process  of  deriving 
a  model  that  best  describes  the  system  dynamics  and  closely  matches  the  actual 
data  presented.  One  method  of  determining  this  model  is  through  least  squares 
estimation.  The  principle  of  least  squares  states  that  the  unknown  parameters  of 
a  model  can  be  deterniined  by  minimizing  a  square  error  criterion.  This  error  is 
made  of  the  sum  of  the  squared  differences  between  the  acutal  signals  and  those 
that  are  computed  or  estimated  [Ref.  1 1 ,  p.  420] .  In  the  generalized  least  squares 
problem,  it  is  assumed  that  the  computed  value  takes  the  form 

y=<|)^(O0 

where  y,<^  are  known  signals  (previous  states/inputs) 

6  are  unknown  parameters 

The  goal  is  to  determine  an  estimate  of  the  system  parameter  matrix  6  such  that 
the  estimated  variable 


y=(0 


(4.4) 
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is  in  close  agreement  with  the  actual  signal,  v.  The  least  squares  method  states 
that  the  parameters  should  be  chosen  in  a  manner  that  minimizes  the  loss  function 
[Ref.  12,  p423] 


W-5) 

4  ,=1 

where  f,  =  v,  -  v,  =  v,  -  -  .  BA 

/  =  1  ,  2,  ...N  N  =  number  of  observations 
The  solution  to  this  problem  is  given  by  the  system  of  equations 

If  the  matrix  0 V  is  nonsingular,  the  minimum  is  unique  and  9  can  be  determined 
by  [Ref.  12,  p.  422], 

For  the  situation  in  which  observations  are  taken  sequentially,  as  is  the 
case  with  a  discrete,  sampled  system,  recursive  equations  may  also  be  derived. 
This  is  known  as  a  recursive  least  squares  algorithm  and  leads  to  recursive 
identification.  This  algorithm  takes  the  previous  observations  and  the  previous 
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least  squares  estimate  for  B  to  calculate  the  current  estimate.  The  least  squares 
estimate  then  satisfies  the  following  set  of  recursive  equations  [Ref.  12,  p.425]. 


z(r)=Jc(0-j:o 

v(r)=M(/)-Mo 


(4.11) 

(4.12) 


where  x(t)  is  the  vector  of  state  variables 

Xq  is  the  set  point 
u(t)  is  the  vector  of  control  inputs 
Uq  is  the  nominal  control  input  at  Xq 
z(t)  and  v(t)  are  the  deviations  from  the  set  points 
Assuming  the  set  point  x^  is  fixed  for  a  particular  operating  condition,  the 
following  relationship  holds. 
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zU)^B  v(i) 


(4.13) 


The  measurable  vector  <p  and  calculated  matrix  6  are  defined  as 
<j)=[Zi(0d:2{0,M{0,::4U),Vi(r),v,(r),v3(r),v4(r),v5(0] 
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The  recursive  least  squares  algorithm  is  utilized  online  to  estimate  the  coefficients 
in  6  from  measurements  of  z(t),  z(t),  and  v(t).  From  this,  A„  and  B„  matrices  trom 
(4.11)  are  derived. 

2.  Linear  Quadratic  Regulator 

Once  we  identify  the  system  parameters  by  least  squares,  many  of  the 
linear,  optimal  control  tools  become  available  for  use  in  control  gain  calculation. 
By  optimal  control  theory  we  can  determine  the  control  input  for  the  purpose  of 
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maximizing  a  measure  of  performance  or  minimizing  a  cost  function.  If  the 
system  we  want  to  control  is  uncertain,  it  is  well  known  that  performance  has  to 
be  sacrificed  for  robustness.  Performance  is  then  measured  by  a  maximization  or 
minimization  of  a  desired  criterion.  This  leads  to  the  concept  of  optimal  stochastic 
control,  a  methodology  which  recognizes  the  random  behavior  of  the  system  and 
attempts  to  optimize  response  or  stability  on  the  average  rather  than  assured 
precision. 

A  discrete  system  is  described  by  the  stale  equation 

The  linear  quadratic  regulator  algorithm  is  designed  to  find  an  optimal  control 
u*(x(k),  k)  that  minimizes  the  performance  measure  [  Ref.  12.  p.  78] 


■k‘0 


(4.16) 


where  0*  is  a  real  symmetric  positive  semi-definite  n  x  n  matrix 

is  a  real  symmetric  positive  definite  m  x  ni  matrix 
The  minimization  of  the  performance  measure  produces  the  control  law  in  the 
form 
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where 


(4.18) 


and  X*  is  an  m-dimensional  adjoint  vector  determined  by 

For  the  indirect  adaptive  control  algorithm,  the  matrix  C*  is  not  a  constant  matrix 
but  is  changing  over  time  in  conformance  with  the  new  parameterization  of  the 
plant. 

The  critical  aspect  of  minimizing  the  cost  function  is  the  proper  choice 
in  the  weighting  of  the  Q  and  R  matrices.  The  larger  the  Q  matrix  in  relation  to 
the  R  matrix  forces  the  minimization  of  the  states  (or  the  variance  from  set  point) 
with  less  emphasis  on  the  control  effort.  If  R  is  proportionally  larger  than  Q,  the 
reverse  is  true.  The  control  utilization  is  minimized  at  the  expense  of  response 
time  and  state  oscillations.  The  desired  control  system  performance  is  thus  dnven 
by  the  selection  of  these  parameters. 


D.  DIRECT  ADAPTIVE  CONTROL 


Direct  adaptive  control,  as  specified  earlier,  combines  the  functions  ot 
parameter  identification  and  control  law  calculation.  For  comparison  and  analysis 
purposes,  the  use  of  a  direct  adaptive  control  algorithm  has  been  investigated.  An 
existing,  generalized  direct  adaptive  control  algorithm,  written  tor  the  estimation 
and  control  of  constrained  multivariable  systems,  has  been  adapted  and 
implemented  for  utilization  in  the  chamber  model.  This  generalized  algorithm 
was  written  by  Dr.  Cory  Finn,  NASA  Ames  Research  Center,  and  can  be  found 
in  Reference  13.  The  following  sections  provide  a  brief  discussion  of  the  direct 
adaptive  control  algorithm. 

1.  System  Identification 

The  state  space  representation  utilized  for  implementation  of  the  direct 
adaptive  control  system  is  referred  to  as  the  autoregressive  moving-average 
(ARMA)  model.  This  is  an  extension  of  the  linear,  discrete  time  state-space  model 


jc(r+l)  =  Ax(t)+Bu(t) 
yit)  =  Cx(f) 


(4.20) 

(4.21) 


which  is  rewritten  as 
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A{q  ‘)v(0  =  B(q'^)u(t-\) 


(4.22) 


where  A(q'‘ )  and  B(q[)  are  polynomials  in  the  backward  shift  operator  q^  .  These 
may  be  expanded  as 


A{q'^)  =  1  +  cl^q'^  +  ....  (4.23) 

B{q-')  =  ^0  +  b^q-^  +  ...  +  bj}'"'  (4.24) 


For  the  chamber  model,  the  parameters  of  the  polynomials  A(q^)  and  B(q  ‘}  are 
unknown.  These  parameters  may  be  estimated  in  terms  of  a  system  model. 

A(,q-^)y{t)=B{q-^)u{l-\) 

The  parameter  estimates  A(q‘)  and  B(q‘)  are  determined  utilizing  a  recursive  least 
squares  algorithm  similar  to  the  one  described  in  the  previous  section.  For 
simplification,  equation  (4.25)  is  written  as 


y{t)=4>^{t)ht) 

where  is  a  matrix  of  the  past  outputs  and  inputs  and  6  is  the  parameter  matrix. 
The  parameter  matrix  is  determined  by  using  a  recursive  least  squares  estimation 
algorithm  with  a  variable  forgetting  factor  [Ref.  11] 


^(04(f-i)+ 


(4.27) 
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Mt)  X(r)^0^(t-\)P(t-\)cp(t-l) 


(4.28) 


e(;)=v(r) l)0fr-l) 


(4.29) 


where  P(t)  is  the  covariance  matrix,  and  eft)  is  the  prediction  error. 
The  algorithm  uses  a  variable  forgetting  factor  which  is  defined  as 


\(/)=X(r-l)  + 


(4.30) 


where  Nq  is  the  memory  length  and  is  a  constant  [  Ref.  14,  p  831-835]. 

2.  Control  Algorithm 

The  generation  of  the  necessary  control  inputs  is  accomplished  via  a 
constrained  multivariable  predictive  controller  [Ref.  13],  The  predicted  output 
signal  at  a  time  step  t+k  can  be  written  as  a  linear  function  of  previous  inputs  and 
output  signals.  Equation  (4-26)  then  becomes 


f(t+k)=a(q'^)y{t)+^{q  ‘iMfr+k-l) 


(4.31) 


where 
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a{q  '')  =dg+d;,  +. 


(4.32) 


(4.33) 

The  polynomials  a(g'')  and  &(q‘)  are  obtained  from  the  polynomials  A(q‘)  and 

^  , 

B(q  )  by  solving  the  Diophantine  equations  [Ref.  11,  p.210]. 

1  =F{q  '^)A{q '')  +q  ‘^Giq  ■‘) 

(4.34) 

F(q-^)  =  l  +  ... 

(4.35) 

G(^'')=8o  ^  8 ^  ••• 

(4.36) 

Given  the  values  of  F(q^)  and  G(q'')  ,  then 

(4.37) 

kq-'>F{q-')6Xq-') 

(4.38) 

It  is  desired  to  solve  the  Diophantine  equations  for  all  values  of  k  between  one  and 
the  prediction  horizon,  T,  The  prediction  horizon  should  be  chosen  such  that  it 
is  larger  than  the  process  dead  time,  which  is  defined  as  the  lag  in  system 
dynamic  response  to  an  input.  This  is  necessary  in  order  to  obtain  stable  control. 
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Th  goal  is  to  solve  the  constained  multivariable  predictive  control  problem  via  an 
objective  function,  minimized  with  respect  to  the  control  input  U 


|K-K‘1  +  wju] 


subject  to  the  constraints 


(4.39) 


u  ^  u  ^  u 

(4.40) 

-6  ^  ^U  ^  6 

(4.41) 

Y^Y^Y 

(4.42) 

where  K*  is  the  vector  of  future  set  points 

t  is  the  vector  of  predicted  future  inputs 
U  is  the  vector  of  future  inputs 

w/  is  the  vector  of  weights  on  the  predicted  future  tracking 
errors 

is  the  vector  of  weights  on  the  input  efforts 
a  is  the  vector  of  lower  bounds  of  the  input  efforts 

V  is  the  vector  of  upper  bounds  of  the  input  efforts 

5  is  the  vector  of  move  size  limitations  on  the  manipulated 
input  variables 

X.  is  the  vector  of  lower  bounds  on  the  output  variables 

Y  is  the  vector  of  upper  bounds  on  the  output  variables 
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A 

For  illustration  purposes,  Y  and  U  are  expanded  as 

(4.43) 

(4.44) 

where  kn  represents  the  number  of  outputs  in  the  process  and  km  represents  the 
number  of  inputs.  For  the  chamber  model,  km  and  kn  are  equal  to  4  and  5. 
respectively. 

The  above  optimization  can  be  cast  as  a  linear  programming  problem, 
and  in  turn  be  solved  using  a  simplex  algorithm  (Appendix  A).  The  set  of  inputs, 
Ui(t),U2(t),...Ukin(t)  calculated  using  the  simplex  algorithm  are  then  implemented 
into  the  chamber  for  control.  The  future  inputs,  u,(t+ 1 u,(t+T-l),... Ukni(t+T- 
1)  that  are  obtained  from  the  optimization  are  never  actually  utilized  for  future 
c.ontrol.  This  is  because  there  is  a  new  set  of  inputs  are  calculated  at  each  time 
t  for  each  time  step. 
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V.  SIMULATION  RESULTS 


A.  STRATEGY 

A  critical  first  test  of  any  preliminary  control  algorithm  is  to  subject  it  to 
rigorous  and  realistic  simulations.  The  model  developed  representing  the  dynamics 
of  the  chamber  air  mass  provides  the  link  for  an  initial  performance  study.  The 
strategy  for  testing  of  the  system  includes  a  simulation  of  the  steady  state 
operational  mode  with  minor  pertubations  in  set  point  and  a  simulated  transition 
or  step  response  from  a  pressure  of  34.464  kPa  (5  psia)  to  52.13  kPa  (7.5  psia). 
The  set  point  requirements  for  various  points  are  provided  in  Table  III.  Finally, 
the  control  system  is  subjected  to  a  seiies  of  incremental  load  increases. 


Table  III  CHAMBER  SET  POINTS 


Pressure 

(kPa) 

O;  Mass 

Fraction 
(Dry  air) 

No  Mass 

Fraction 
(Dry  air) 

Temperature 

(“C) 

Absolute 

Humidity 

(kg  H;0/kg  dry  air) 

101.325 

0.231574 

0.768426 

21.1 

0.0079504 

86.16 

0.272867 

0.727133 

21,1 

0.0093431 

68.92 

0.338166 

0.661834 

21.1 

0.0116668 

51.69 

0.444653 

0.555347 

21.1 

0.015537 

34.46 

0.64903 

0.35097 

21.1 

0.0232739 
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1.  Steady  State  Response 

In  this  section,  we  addressed  the  performance  of  only  the  author  derived 
indirect  adaptive  controller/algorithm.  A  system  simulation  representing  nominal 
steady  state  loads  has  been  conducted.  This  loading  is  defined  in  Table  II.  The 
nominal  control  inputs  necessary  to  maintain  the  set  points  were  derived  utilizing 
mass  and  energy  balances  over  the  entire  system.  Appendix  B  contains  the  results 
of  the  parametric  studies  utilized  to  obtain  steady  state  values  for  the  condensing 
heat  exchanger  performance.  Simple  mass  balances  were  used  to  determine  the 
necessary  oxygen  and  nitrogen  makeup  flows  in  addition  to  the  required  vacuum 
flow. 

The  simulation  is  initiated  by  varying  the  control  inputs  randomly  within 
a  10  percent  range  of  the  steady  state  values.  This  was  accomplished  for  two 
reasons;  (1)  to  allow  the  recursive  least  squares  algorithm  to  "learn"  the  system 

A 

model  and  allow  the  estimated  parameter  matrix  6  to  converge,  and  (2)  to  allow 
the  system  to  drift  from  the  initial  set  point.  Three  hundred  time  steps,  equivalent 
to  1.25  hours  of  simulation,  has  been  alloted  for  the  converging  period.  At  this 
juncture,  the  controller  was  turned  on.  The  dynamic  response  of  the  system  is 
represented  in  Figures  5.1,  5.2,  5.3  and  5.4.  When  the  control  system  is  "turned 
on",  the  states  fluctuate  slightly  but  remain  within  the  limits  of  control.  The 
fluctuations  that  can  be  observed  for  the  first  hour  in  the  plot  in  Figure  5.1  is 
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representative  of  the  system  "learning”.  The  least  squares  estimator  managed  to 
recursively  converge  relatively  quickly  and  presents  an  excellent  model  of  the  plant 
at  this  operating  condition.  The  controller  response  to  the  minoi  set  point 
pertubation  and  adaptation  to  the  steady  state  scenario  is  satisfactory.  However, 
a  steady  state  error  is  observed  in  all  states.  The  largest  error  encountered  is  in 
temperature.  The  temperature  settled  into  a  steady  state  value  of  2[A°C,  which 
translates  into  a  steady  state  error  of  1.4%. 


Figure  5.1  Pressure  Steady  State  Response  at  34.464  kPa 
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Figure  5.2  Oxygen  Mass  Fraction  Steady  State  Response 


Figure  5.3  Temperature  Steady  state  Response 
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Chwnber  Homiditjr 


Figure  5.4  Absolute  Humidity  Steady  State  Response 

2.  Step  Response 

In  addition  to  stability,  sensitivity  and  accuracy,  we  are  always 
concerned  with  the  transient  response  of  a  feedback  system.  Transient 
characteristics  are  normally  defined  on  the  basis  of  a  step  response.  Therefore, 
the  next  performance  evaluation  undertaken  is  a  simulation  of  the  system  in  a 
transition  from  a  state  defined  at  34.464  kPa  to  one  defined  at  51.69  kPa. 
Transition  control  in  a  hypobaric  environment  is  critical.  The  key  to  this  transition 
lies  in  the  ability  of  the  system  to  maintain  the  mass  fraction  (or  partial  pressure) 
of  oxygen  in  the  proper  proportion  to  the  total  mass  (or  total  pressure).  The 
physiological  effects  of  an  improper  transition  to  lower  pressures  are  great  and 
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must  be  addressed  in  any  control  algorithm.  For  the  transient  response,  both  the 
direct  and  indirect  adaptive  controllers  will  be  investigated. 

a.  Indirect  Controller  Performance 

The  adaptive  based  linear  quadratic  regulator  controller  performance 
is  provided  in  Figures  5.5,  5.6,  5.7  and  5.8.  The  pressure  and  oxygen  mass 
fraction  transitions  are  adequate  (Figures  5.5  and  5.6).  The  transitions  are  both 
essentially  linear  and  achieve  the  new  steady  state  at  approximately  the  same  time 
with  minimum  over  shoot.  This  is  an  indication  that  the  proper  proportion  of 
oxygen  is  maintained  during  the  transition.  Since  an  upward  transition  in  pressure 
does  not  involve  any  physiological  considerations  (i.e.,  the  bends)  and  there  are 
no  current  experimental  time  constraints  imposed,  the  settling  time  is  not 
considered  a  critical  performance  measure.  The  main  deterrent  to  more  rapid 
transition  is  the  physical  constraints  imposed  by  the  makeup  flow  rates.  The 
temperature  transient  response  is  presented  in  Figure  5.7.  A  comparatively  large 
jump  is  observed  at  the  onset  of  the  set  point  change.  This  may  be  due  in  part  to 
the  coupling  of  the  state  variables.  Pressure  and  oxygen  mass  fraction  carry  a 
much  greater  weight  in  the  control  law  development,  ie.  the  Q  matrix  in  the  LQR 
algorithm  is  weighted  significantly  towards  these  variables.  As  a  result,  the 
controller  initiated  an  immediate  increase  in  pressure.  In  conformance  with  the 
ideal  gas  law  for  a  fixed  volume,  an  increase  in  pressure  is  accompanied  by  an 
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increase  in  temperature.  Sines  the  temperature  set  point  did  not  change  in  the 
transition,  its  initial  response  is  only  a  result  of  the  changes  initiated  for  humidity 
and  pressure.  The  humidity  transition  to  the  new  set  point  is  adequate  with  a 
minor  spike  encountered  at  the  outset,  largely  due  to  an  initial  jump  in 
temperature.  Numerous  follow  on  simulations  were  conducted  to  reduce  this 
temperature  anomaly  by  adjusting  the  weighting  of  the  Q  and  R  matrices  with  little 
success.  The  steady  state  errors  at  the  new  set  point  exhibit  minimal  variance 
from  the  errors  observed  at  the  initial  set  point. 


Figure  5.5  Pressure  Transient  Response 
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b.  Direct  Controller  Performance 

The  direct  control  algorithms  are  written  in  Fortran  and  the 
computer  code  is  presented  in  Appendix  A.  These  have  been  integrated  into  the 
Matlab  model  for  simulation  and  analyses.  The  transient  response  of  the 
predictive,  multivariable  control  algorithm  to  a  change  in  set  point  is  presented  in 
Figures  5.9,  5.10,  5.1 1,  and  5.12.  The  responses  are  very  similar  to  those 
obtained  from  the  LQR  controller.  The  transitions  in  pressure  and  oxygen  is 
almost  identical  in  both  shape  and  settling  time  (Figures  5.9  and  5.10).  There  was 
slightly  less  overshoot  in  pressure  transient.  In  terms  of  temperature,  the 
pertubations  about  the  set  point  are  comparatively  'mailer  (Figure  5.11).  A  spike 
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of  only  about  I  “C  was  observed,  compared  to  over  2  “C  for  the  indirect 
controller.  The  humidity  response  in  the  direct  algorithm  represents  an 
improvement  over  the  indirect  response,  largely  due  to  the  minimization  of  the 
temperature  variance. 

The  steady  state  errors  for  the  new  set  point  are  near  zero.  This 
can  be  accounted  for  in  part  to  the  computationally  intensive  simplex  and 
Diophantine  algorithms.  The  30  hour  simulation  required  over  24  hours  of  actual 
computer  time  on  a  VAX  system,  compared  to  only  about  15  minutes  for  the  LQR 
based  controller  on  a  IBM  80486  personal  computer.  The  computational  intensity 
of  the  direct  controller  severly  limits  it  viability.  Computer  code  streamlining  and 
optimization  would  be  required  for  implementation. 


65 


ftaatawk(  *7  •* 


aTO 

0^ 

<M0 

0J3 
I  OJO 

I 

0.4S 

040 

I 

Figure  5.10  Oxygen  Mass  Fraction  Transient  Response 


Figure  5.11  Temperature  Transient  Response 
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Figure  5.12  Absolute  Humidity  Transient  Response 

3.  Load  Disturbance 

An  effective  manner  in  which  to  test  the  adequacy  of  the  control 
algorithm  is  to  subject  it  to  demanding  loading  conditions.  Since  the  range  of 
potential  loads  in  the  CERC  can  vary  significantly  depending  upon  the  prevalent 
human  acitivity,  it  is  necessary  to  subject  the  controller  to  a  large  variation  in 
loading.  Hence,  the  response  of  the  LQR  based  controller  to  a  series  of  step 
disturbance  increases  was  analyzed.  For  this  analysis,  the  latent  and  sensible  heat 
loads  in  additon  to  the  oxygen  consumption  have  been  increased  incrementally  over 
a  50  hour  simulation  at  the  same  operating  set  points.  A  graphical  representation 
of  the  variable  loading  is  available  in  Figures  5.13,  5.14  and  5.15. 
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In  order  to  accurately  appraise  the  performance  of  the  closed  loop 
controller,  it  is  necessary  to  analyze  the  open  loop  system  response  to  the 
incremental  disturbance.  This  response  is  presented  in  Figures  5.16,  5.17,  5.18 
and  5.19,  The  nonlinearity  of  the  dynamic  system  can  be  readily  observed  in  the 
first  order  response. 
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5.16  Pressure  Open  Loop  Response 


Figure  5.17  Oxygen  Mass  Fraction  Open  Loop  Response 
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5.18  Temperature  Open  oop  Response 


Figure  5.19  Absolute  Humidity  Open  Loop  Response 
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The  load  disturbance  on  the  closed  loop  system  is  depicted  in  Figures 
5.20,  5.21,  5.23  and  5.24.  Looking  specifically  at  the  pressure  response  (Figure 
5.20),  it  is  observed  that  only  a  series  of  spikes  corresponding  to  the  onset  of  the 
disturbances  deters  from  a  relatively  flat  curve.  The  maximum  pertubation  from 
set  point  only  represents  a  1.4  percent  variance.  In  light  of  the  magnitude  of  the 
load  disturbance  from  assumed  steady  state,  the  response  may  be  classified  as 
adequate. 

The  oxygen  mass  fraction  response  is  excellent  (Figure  5.21).  Despite 
tripling  the  oxygen  consumption  and  latent  heat  production,  the  largest  observed 
deviation  from  set  point  is  just  0.3  percent.  Temperature  variation,  depicted  in 
Figure  5.22,  is  comparatively  large.  The  relative  weighting  of  the  temperature, 
as  discussed  previously,  is  one  probable  cause  of  the  larger  deviation  from  set 
point.  Efforts  to  improve  the  temperature  response  by  adjusting  the  weight  had 
a  tendency  to  cause  the  controller  to  cycle  and  go  unstable.  Despite  over  doubling 
the  overall  latent  and  sensible  loads,  the  maximum  variation  was  still  limited  to  2 
"C.  The  tracking  and  control  of  humidity,  because  of  the  strong  coupling  with 
temperature,  had  a  similar  response.  The  maximum  overshoot  was  observed  to 
be  0.002  kg  water/kg  ory  air,  or  about  eight  percent. 
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Figure  5.21  oxygen  Mass  Fraction  Disturbance  Response 
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Figura  5.22  Temperature  Closed  Loop  Response  to  Disturbance 


Figure  5.23 
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VI.  CONCLUSIONS/RECOMMENDATIONS 


A.  SUIVIMARY 

A  preliminary  control  algorithm  and  mathematical  model  have  been 
developed  for  the  atmospheric  system  within  the  Controlled  Environment  Research 
Chamber  located  at  NASA  Ames  Research  Center.  The  basic  system 
configuration  and  groundwork  parametric  studies  have  been  reviewed  for  inclusion 
in  the  basic  modeling  equations.  Adaptive  control  techniques  have  been  developed 
and  tested  for  potential  utilization  in  system  identification  and  control  of  this 
nonlinear,  multivariable  process.  Included  within  these  techniques  are  algorithms 
involving  recursive  least  squares  estimator  for  the  approximation  of  system 
dynamics. 

The  indirect  adaptive  control  scheme,  utilizing  a  linear  quadratic  regulator 
for  control  law  formulation,  provided  satisfactory  control  of  the  chamber 
atmosphere.  The  controller  effectively  maintained  system  state  variables  in  a 
steady  state  simulation  with  minimal  steads  state  error.  Additionally,  the 
algorithm  tracked  a  set  point  change  with  essentially  no  overshoot  and  a  finite 
steady  state  error.  The  introduction  of  significant  metabolic  step  disturbances  in 
system  operation  caused  only  small  deviations  in  the  atmospheric  set  points. 
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The  direct  control  algorithm,  adapted  in  cooperation  with  Dr.  Cory  Finn, 
exhibited  a  greater  degree  of  accuracy  in  set  point  maintenance.  The  transient 
response  to  a  set  point  change  was  similar  in  form  to  the  indirect  algorithm. 
However,  it  was  computationally  intensive  and  in  us  present  configuration  is  not 
a  viable  alternative. 

Presuming  the  model  developed  adequately  portrays  the  chamber  dynamics, 
adaptive  control  techniques  appear  to  be  an  excellent  alternative  for  the 
simultaneous  control  of  temperature,  pressure,  composition  and  humidity. 


B.  RECOMMENDATIONS 

Significant  additional  work  is  necessary  prior  to  implementation  of  this 
adaptive  control  algorithm.  The  following  is  a  listing  of  areas  which  will  require 
attention: 

•  The  mathematical  model  developed  requires  the  inclusion  of  the  dynamics  of 
the  control  system  actuators,  measurement  and  metering  devices  and 
associated  time  delays.  In  addition,  mass  and  energy  diffusion  dynamics 
should  be  analyzed  for  potential  inclusion  into  the  model. 

•  Strawman  experiments  involving  chamber  utilization  must  be  refined  to 
accurately  define  control  system  requirements.  These  experiments  would 
provide  transient  response  characteristics  required  as  well  as  a  better 
representation  of  overall  operating  envelopes. 

•  Increased  work  on  hardware  selection  and  parametric  performance  studies. 
In  particular,  a  suitable  condensing  heat  exchanger  must  be  selected  and 
performance  criteria  implemented  into  the  chamber  dynamic  model. 
Additionally,  performance  data  needs  to  be  refined  for  the  gas  supply  system. 
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•  Additional  simulations  need  to  be  O'-'nducted  on  the  improved  model  to 
optimize  performance  of  the  algorithm.  Eventual  implementation  into  the 
chamber  requires  exhaustive  utilization  of  preliminary  testing  to  maximize  the 
performance. 
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APPENDIX  A.  COMPUTER  PROGRAMS 


I .  This  first  section  details  the  model  and  the  indirect  adaptive  control  routine 
that  utilizes  an  LQR  control  gain  calculation. 


%  Program  cerc  .m 


lie!  bhml  .met 
!dei  bhmi.met 
'del  bhmd.met 
'del  bhm4.met 
fig; 
clear: 

. . . . . . 

%  This  computer  code  represents  a  dynamic  model  and  indirect  adaptive  control 
%  rouiineot'the  temperature,  composition,  pressure  and  humidity  of  the  habitat 
%  enclosedwithin  the  Closed  Environment  Research  Chamber  that  is  currently 
%  undergoingrenovationat  NASA  Ames  Research  Center  located  at  Moffet  Field, 
%  CA 

**«••*««««*»««««•«««•«««»•*•««*«•««««»«««*«««•••*•««■**»• 

%  The  variables  utilized  for  state  definition  and  control  are  as  follows: 

"c  xl  =  Chamber  pressure  (kPa) 

X.2  =  Mass  Fraction  of  02  in  dry  air 
%  x3=  Chamber  Temperature  (degrees  C) 

%  x4=  Chamber  absolute  humidity  (kg  H20/kg  dry  air) 

%  ul  =  Vacuum  flow  rate  out  of  chamber  (kg/hr) 

%  u2=  Mass  flow  rate  of  oxygen  into  chamber  (kg/hr) 

%  u3=  Mass  flow  rate  of  nitrogen  into  chamber  (kg/hr) 

%  u4=  Coolant  stream  entrance  temperature  into  CHX  (degrees  C) 

%  u5=:  Air  flow  rate  from  chamber  into  CHX  (kg/hr) 

%  The  following  are  variables  that  are  utilized  in  system  defintion: 

%  U  =  Heat  exhanger  performance  parameter  iKI/m*2-C) 

%  A=  Heat  exchanger  area  (m*2) 

%  mc=  Coolant  flow  rate  across  CHX  (kg/hr) 

%  R=  Universal  gas  constant  8.314  kj/kgmole  K 
%  Vc=  Chamber  volume  (m‘3) 

%  mieak  =  Leakage  of  ambient  air  into  the  chamber  (kg/hr) 

%  xl=  Mass  fraction  of  02  in  ambient  (outside)  air 
%  02  =  Human  consumption  of  02  (four  persons)  (kg/hr) 
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%  H20=  Human  produciion  ot  H20  (tour  personsMkg'hri 
't  qh=  Human  sensible  heat  load  (four  personukJ  hn 

^m=  Machinerv  miscellaneous  sensible  heat  loads  lU  hn 
'c  a  1  .a2.a3 ,b  1  .h2.b3  =  Coefficients  ot  heat  capacities  tor 
chamber  constituents  iO2.N2,H20l 

'■f  vap=  heat  ot'  vaporization  ot  H20  from  human  m  cham'oer-2287  U  hr 

*'{  hvap  I  =  Heat  of  vaponzation  of  H20  in  CHX  exit  stream 

%  0=  Matnx  consunt 

‘t  \3  =  Chamber  temperature  in  degrees  K 

‘i  x4=  CHX  air  exit  temperature  in  degrees  K. 

"  The  simulation  of  the  multivariable,  non  linear  equations  will  be 
set  up  as  follows: 

xdot=  Fixl’xdot  -f  G(x.u)  -e  Hlx.w); 

^  Solving  for  xdot  , 

"c  xdot=  iI-F<x)l-l*G(x,U)  -r  \\  -  FixD-l*  Hix.wi 

"c  where  x  is  a  vector  of  state  vanables.  u  is  a  vector  of  the  control  inputs 
%  and  w  is  a  vector  of  the  metabolic  loads. 

%  Setting  the  initial  conditions  of  the  states  and  control  variables  for  a  simulation 
%  involving  a  state  defined  at  34.464  kPa.  In  addition,  the  set  point  at  this  sute 
%  is  given  as  well  as  the  maximum  control  inputs. 

x=:134.464;0.64903;21.1;0.023273911; 
xs«t  =  |34.464;0.64903;21 . 1 .0.0232739); 
u  =  ll.0',l.0;0.3658.3;980); 
uset  =  ll;I.00;0.3658;1.5;9801; 
umax  =  110;10;10;4;2000|; 

%  Setting  the  other  variables 

%  Estimated  volume  of  the  chamber  (m3) 

Vc=157; 

%  Heat  capacity  coefficients  (1-02.  2-N2.  3-H20) 

a  I  =0.9324; 

a2  =  0.90468; 

a3=l  83T. 

hi  =0.000184; 

b2  =  0.000204; 

b3  =  0.000744; 

%  Machinery/miscelianeous  sensible  heat  load  inside  chamber  (kJ/hr) 
qm  =  8586.33; 

%  Oxygen  consumption  of  humans  (4  persons-3(@normal  activity, 

%  i(®exercis«)  (kg/hr) 

02=0.3956; 

%  Estimated  Latent  heat  production  ■  humans  i  water)  (kg/hri 
H20=  0  7398; 

Estimated  Sensible  heat  production  -  humans  (U/hr) 
qh  =  2213; 

%  Standard  heat  of  vaporization  of  water  inside  chamber  i@  21 .  loC) 

■•ap  =  2287; 
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^  Estimated  air  leikaae  into  chamber  ikg/hr) 
mleak =0.0025: 

Oxygen  mass  traction  ot  air  leaking  into  the  chamber 
xi=0  2344; 

‘'r  Universal  gas  constant  ikJ/kgmole-Ki 
R  =  8  314; 

‘Z  Initializing  the  vectors  for  utilization  in  Uie  recursive  least  squares  algonthm 

alpha=0.0001; 

1)1  =alpha*eyei9); 
thl  =zeros(9.1); 
q2  =  alpha*eye(9); 
th2= zeros!  9,1); 
q3  =  alpha*eye(9); 
th3=zeros(9.1); 
q4  =  alpha*eye(9); 
th4=zerost9.1); 


7^,,**. a...*.*..*....*. 

%  Setting  the  overall  simulation  length,  integration  time  step,  and  sampling  time 

%  Overall  number  of  sampling  time  steps 
n=2000; 

%  Sampling  instance  (i.e.  p  integration  time  steps  per  sample) 

p  =  6-. 

%  Integration  lime  step  (hrs) 

%  T  is  the  sampling  duration  required  in  the  integration  of  the  equations 
%  of  motion  that  define  the  dynamics  of  the  processes  involved  in 
%  regulating  temperature,  humidity,  pressure  and  composition- 
%  Simulation  for  n  integration  steps 
T=0.004I67; 

%  Setting  up  the  time  vector 
time  =  (l;n*p+  l|; 
time=time*T; 

%  Running  the  simulation  for  n  samples 
for  k=  I  :n; 


%  The  simulation  will  run  for  p  integration  time  steps  prior  to  sampling 
%  and  evaluating  the  data  for  parameter  estimation  and  control  input. 

p  =  6; 

r  =  p-l; 

for  i  =  p*k-r;p*k 


xl  =x(l.i); 
x2  =  x(2,l); 
x3  =  x(3,i); 
x4=x(4.i); 
x5  =  x(3.i)+273; 
xe  =  u(4, 0  +  281.33; 
x7  =  u(4.0  +  8.33-. 
ul  =u(l,i); 
u2  =  u(2,i); 
u3  =  u(3,i); 
u4  =  ut4,i); 
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u5  =  u(5.i); 


%  Change  in  the  set  points  at  t  =  5  hours 
%  if  i>  1200 

%  Nset  =  1 5 1  .69;0,444653  ;2 1 . !  ;0  .0 1 553  7); 

%  end; 


Incremental  increases  in  metabolic  loads 

if  i  >  3000 
%  H20=  I : 

%  qh  =  3000; 

%  02=0.8; 

%  end; 

%  if  i  >6000 
%  H20=l.5. 

%  qh=4000; 

%  02=1.1; 

%  end; 

%  if  1  >9000 
%  H20=i.7; 

%  qh  =  5000; 

%  02=1.5; 

%  end: 

^ «*»*«««***»*  ■*««•«**«««*«««*««««**««««*««««««*«•«•«***«****««*«*«•« 

%  Calculate  the  partial  pressure  of  H20  in  the  exit  stream  of  the 
%  CHX 

Psur=  10*(28.5905-8.2/2.306»log(x6)  +  0.002484»x6-  3142/x6); 

if  x3  <x7 

Pstar=0.0; 

end; 

%  Calculate  the  heat  of  vaporization  of  the  moisture  in  the  exit 
%  stream  of  the  CHX 

Hv  =  2213; 

hvapl  =|13.1 1  -  log(Pstar*0.9869f)*Ryi8.02*x6; 

%  Calculate  the  hunudity  of  the  CHX  exit  stream 

we=  18.02/(4*x2-r28+  18.02*x4)*(100*PsUr/(xl-100»PsUr)); 


%  Calculate  me  enthalpies  of  the  entrance  and  exit  gases  of  the  CHX 

CpTl  =x2*{al*x3  +  bl»x3*2)  +  (I-x21*(a2*x3  +  b2*x3*2)  ...; 

+  x4*t  a3*x3  +  b3»x3*2); 

CpT2  =  x2*(al*x7  +  bl»x7*2)  +  (l-x2)*(a2*x7  +  b2»x7*2) .... 

+  we»(a3»x7  =  b3»x7*2); 

%  Calculate  the  change  in  sensible  heat  removed  across  the  CHX 
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Jelh  =  CpT2-  CpTl; 


%  Calculate  the  chanae  in  absolute  humidity  across  the  CHX 

delw  =  we  -x4; 

It  delw  >0 
Llelw=0.0: 

end; 

%  Calculate  the  latent  heat  removed  across  the  CHx 

dell=delw»hvapl; 
delh20  =  u5/(l  -rx4)*del'w. 

%  Calculate  the  total  heat  removed  across  the  heat  exchanger 
delchx  =  u5/(  1  -r  x4)*(delh  +  dell); 

%  Calculate  the  function  Q 

Q  =  (x2*(al  -  al*x3/x5  +  2*hl*x3  •  bl*x3''2;x5)  *  tl-x2)*ia2  -  ; 

a2»x3/x5  +  2*b2*x3  -  b2*x3‘2/x5>  -r  x4«(-  Hv/x5  +  a3  -  a3*x3/x5  * 
2*b3*x3  -  b3*x3*2/x5)); 


%  The  followtng  code  calculates  the  change  in  the  state  variables  given  previous 
%  states  and  current  input 

%  The  following  is  the  mathematical  representations  of  the  matrices 
%  F,  G.  and  H 


F=|0,0.xl/x5.0;  -x2/xl.0,x2/x5.x2/(I +x4);-(x2*(al* 

x3  +  bl*x3*2)  +  (l-x2)*(a2*x3  +  b2*x3"2)  -l-  x4»(Hv  +  a3*x3  +  .... 

b3*x3*2))/(xl*QJ.-(al*x3  +bl»x3*2  -  a2*x3  -  b2»x3*2V  ...; 

Q.  0,-((-i/fl  +  x4))«(x2*(al»x3  +  bl*x3'2)+  <1  -  x2)*  .  .; 

(a2*x3  +  b2*x3'2))  +  Hv  +  a3*x3  +  b3»x3*2  -  (x4/(l  +  x4))»(Hv  +  a3*...; 
x3  +  b3*x3*2))/Q; -l/(xlM/x4-  l/(x4+ 1))),  0,  l/(x5*  ...; 
fl/x4  -l/fx4+l))).0|; 

G  =  (R*x5/<{4*x2  +  28  +  x4*l8)»Vc)*{u2  +  u3  -  ul  +delh20  ); 

(1  +  x4)*R*x5/((4*x2  +  28  +  l8*x4)»Vc*xl)*(u2-  x2*  ...; 
ul);  (I  +x4)*R*x5/{(4»x2  +  28  +  18»x4)*Vc*xl*Q)*delchx;  ..; 

(1  +  x4)*R«x5/((4*x2  +  28  +  18.02*x4)»xl»x4*Vc»(l/x4..  ; 

1-  l/(x4+  l)))*u5/(l +x4)*delw|; 

H  =  fR*x5/((4»x2  +28+  18»x4)*Vc)»(-02  +  H20  +  mleak);  R*x5»(l  +  x4)  ; 

/((4»x2  +28+  18*x4)*Vc*xl)»(xl*tnleak-  02);  R*(l  +  x4)»x5/((  ..; 

4»x2  +28+18*x4)*Vc*xl*Q)»(qh+  qm  +  H20*vap>;  R*x5*(l  +  x4)/((  .... 
4»x2  +28  +18*x4)*Vc*xl*x4»(l/x4-  l/{x4  +  1)))«H20I; 

%  The  following  sequence  will  calculate  the  change  in  x  over  the 
%  integration  time  period 

xdot  =  inv(eye(4)-F)*(G  +  H); 
x(:.i+ l)  =  x(:.i)  +  xdot*T; 
u(:.i+l)=u(:,l); 
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end 

"c  The  t'ollowing  code  will  calculate  the  variables  utilized  m  the  recursive  least 
squares  algorithm  and  will  vary  the  steady  state  control  inputs  tor  the  t'irst  300 
Integration  time  steps  in  order  for  the  parameter  matrix  to  converge 

It'  i<300 
z  =  xi:,i-r  l)-xset; 

V  =  ui :  .il-uset; 
phi  =  lz'.vT; 
zdot'  :.i)  =  xdof. 

%  The  following  code  calls  the  recursive  least  squares  subroutine  for 
|th  1  .q  i  I  =  rlslth  1  .q  I  .phi.zdotl  1  .i)): 

(ih;.q2|  =  rls(th2,q2.phi.zdot(2.i)); 
lth3.q3|  =  ris(lh3.q3.phi.zdot(3.i)); 

[th4,q4|  =  rls(th4,q4.phi,zdot(4,i)); 

wir.l)  represents  the  diffference  between  the  estimated  states  and  the  actual 
%  output  signals 

v  I  :,i)  =  iphi'*|th  l.th2.th3.th4|)'-xdot; 

^  Random  variation  in  steady  state  inputs 

ui  i.i-r  l)  =  u(l>+0.05*rand; 
u(2,i+  l)=u(2)+0.05*rand; 
u(3,i  +  l)  =  u(3)+0.05*rand; 
u{4,i+  l)  =  u(4)+0.01*rand; 
u(5,i+  n=u(5)  +  5*rand; 

end 

%  The  controller  is  turned  on  at  integration  time  step  greater  than  300,  The 
%  recursive  least  squares  parameter  estimator  is  left  on  line  for  contimued 
%  system  identification. 

if  i  >  300 
z  =  x(;,i-r  l)-xset; 
v  =  u(:,i)-uset: 
phi  =  |z’.V|’; 
zdot  =  xdot; 

f  th  1  .q  1 1  =  risdh  I  .q  1  .phi.zdotl  1 )); 
lth2  .q21  =  rls(th2.q2.phi.zdot(2)): 

(th3  .q3 1  =  rls<th3  .q3.phi,zdot(3)); 

( th4.q4|  =  rls(th4,q4.phi.Zdot(4)); 
a  =  |thl{l:4)-; 
th2n;4)’; 
th3(l:4T; 
th4(l:4)'l; 
h  =  lthl(5:9)'; 
th2(5:9)  ; 
th3(5:9)’; 
th4(5:9i’l-. 

%  Linear  quadratic  Regulator  implementation 

%  The  q  vector  is  for  the  sutes  relative  weighting  and  c  vector  for  the  input 
%  relative  weighting 

q  =  3  ‘(T.O.O.OtO,  1 800.0.0;0,0, 10.0;0, 0,0,40); 

c  =  1 .0 1 ,0,0,0,0;0,  .0 1 ,0,0,0;0,0. .  1 ,0,0;0,0, 0,0.0 1 .0;0.0.0,0,0. 10); 
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m  =  lqr(a,b.i),ci; 

^  Control  law  implementation  where  m  is  the  cam  matnx 
iii:.i-r  l)  =  -m*z; 


111 :  .1  1 1  ui . .  1 1  r  ui  :.i  —  1 ), 

%  Constraints  on  calculated  inputs 

if  u(l,i-rll<0 
ui!.i+  l'i=0.0, 
end 

if  u(2.i-r  1)  <  0 
u(2.if-i)=0.0; 
end 

if  u(3.i  f  1 1  <0 
iit3.i  i )  ==0.0; 
end 

if  u(4.i+  ll<0 
u(4,i  +  l)  =  0.0: 
end 

if  u(5.i+  1)  <0 
u(5,i-t- 1)=0.0; 
end 

ifu(l.i+l)>umax(l) 
u{l  ,i  +  l)=sumax(l); 
end 

if  u(2.i+  h>umaxi2) 
u(2,i->-  l)  =  umax(2); 
end 

if  u(3.i+  1)  >umax(3) 
u(3.i+  l)  =  umax(3); 
end 

if  u(4.i+  I )  >iimax(4) 
ut4.i+  i)  =  umax<4); 
end 

if  u(5,i  -r  I )  >umax(5) 
u(5.i+  1  )  =  umax(5); 
end 

end; 

end; 


%  Plots  of  states  versus  time 

ploUtime.xfl.:)); 
xlabeKTime  (hours)’) 
ylabelC Pressure  (IcPa)') 
titleCChamber  Pressure') 
meta  bhtn  I ; 
pause; 


84 


plot(time.xt2.:)), 

xIabeU'Time  (hours)') 

ylabiili''  ^  .-.s  Fraction  ot' 02  iDry  Airi'i 

titlet’Mass  Fraction  ot  02') 

meta  hhm2; 

pause; 

plot(time.x(3.;ii, 
xlabeli'Time  (hours)') 
vlabeirTemperature  (O’) 
title  1 'Chamber  Temperature' 1 
meta  bhm3; 
pause; 


plot(time,x(4.:)); 
xlabeli  'Time  (hours)') 

> label!' Absolute  Humidity  ikg  H20i'ktf  dr>  ain't 
title) 'Chamber  Humidity  ) 
meta  hhm4. 
pause; 


%  Plots  of  inputs  versus  time 

xiabelCTime  (hours)’) 
ylabeU’Vacuum  Flow  Rate  (kg/hr)') 
title! 'Vacuum  Flow  Rale’) 
pause; 

plot(time.u!2.;)) 

xiabelCTime  Ihoursi') 

ylabel! 'Oxygen  Flow  Rate  (kg/hr)') 

title! 'Oxygen  Flow  Rate') 

pause: 

plot(llme,u(3,:)) 

XiabelCTime  ihours)') 
ylabelCNitrogen  Flow  Rate  (kg/hr)') 
tilleCNitrogen  Flow  Rate’) 
pause. 


plo!!!ime.u(4. ;)) 
xiabelCTime  ihoursi’) 
ylabel! 'Temperature  iC)') 
titlel'CHX  Entrance  Temperature’) 
pause: 

plotItime.uiS,:)) 
xiabelCTime  (hours)’) 
ylabelCAir  Flow  !kg/hr)') 
titleCCHX  Air  Flow  Rate’) 
pause: 
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**************>|c**>,.******:«***********;t:**>K**!*:*:t:****:^**:t.****:t:***#:=t=*** 

tiinction  Ithnew,  Onew|  =  ris(thoid,  Ot>ld.  phi.  zdoil 

%  tunction  Ithnew,  Onewt  =  rlsithold.  Qold,  phi.  zdoti 

"c  Recursive  Least  Squares  for  estimating  the  rows  ot  the  A  and  B  matrices. 

"  th  =  parameter  vector  (in  our  case  9  by  1 ) 

^  0  =  positive  definite  matnx  (9  by  9i  Initialize  ii  at  0  =  aipha*eyei9i. 
with  alpha  a  small  positive  number  isey  0.0001 1 

"c  phi  =  9hy  1  vector  =  (zld) . z4(t.).vl(ti . v5(or 

Qnew  =  Qold  +  phi*phi ' ; 

P=inv(Qnewi; 

thnew  =  thold  -h  P*phi*(zdoi-phi’»thold) ; 
end  %  rls 

*5le*:(£**HcHr**:(r*:(::ts****;i;***********Hcit:****:t:**:*::*:***5(c**!(e5t:**:(e**5t:*********** 

2.  The  second  section  presents  the  model  and  Fortran  computer  code  utilized 
in  the  development  of  the  direct  adaptive  control  routine  [Ref.  13], 

a.  Modified  Matlab  Model 

^  Program  cerc.m 

%  This  program  is  a  simulation  of  the  temperature,  composition.pressureand 
%  hurnidityof  a  Controlled  Environment  Research  Chamber  that  is  curremly 
%  undergoing  renovation  at  NASA  Ames  Research  Center  located  at  Moffet  Field 
%  CA 


%  The  variables  utilized  for  sute  definition  and  control  are  as  follows; 

%  x(  l)  =  Chamber  pressure  I LPa) 

%  x(2)  =  Mass  fraction  of  02  in  dry  air 

%  x(3)  =  Chamber  temperature  (degrees  C) 

%  x{4)  =  Chamber  absolute  humidity  (kg  H20/kg  dry  air) 

%  u(l)  =  Vacuum  flow  rate  out  of  chamber  (kg/hr) 

u(2)  =  Mass  flow  rate  of  oxygen  into  chamber  (kg/hri 
%  u(3)  =  Mass  flow  rate  of  nitrogen  into  chamber  ikg/hri 

%  u(4)  =  Coolant  stream  entrance  lemperacure  into  CHX  (degrees  C) 

%  u(5)  =  Air  flow  rate  from  chamber  into  CHX  (kg/hr) 

%  The  following  are  variables  that  are  utilized  in  system  definiion; 

%  U  =  Heat  exhanger  performance  parameter  (KJ/m*2-C) 

^  A  =  Heat  exchanger  area  (m'2) 

%  me  =  Coolant  flow  rate  across  CHX  (kg/hr) 

%  R  =  Universal  gas  constont  8.314  kJ/kgmole  K 

%  Vc  =  Chamber  volume  (m'3) 

^  mleak  =  Leakage  of  ambient  air  ituo  the  chamber  (kg/ht) 

%  xl  =  Mass  fraction  of  02  in  ambient  (outside)  air 
%  02  =  Human  consumption  of  02  (four  persons)  (kg, hr) 

%  H20  =  Human  production  of  H20  (four  persons)  (kg, hr) 

%  qh  =  Human  sensible  heat  load  (four  person)  (kJ/hri 

%  qm  =  Machinery/miscellaneous  sensible  heat  loads  (kJ, 'hr) 

%  vap  =  Heat  of  vaponzation  of  H20  from  human  in  chamber-2287  U/hr 
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"c  hvapl  =  Heat  of  vaporization  ot  H20  in  CHX  exit  stream 
%  hv  =  Heat  ot  vaporization  ot' H20  - 2213  kJ'hr 

"r  Q  ~  Matrix  constant 
"  x5  =  Chamber  temperature  in  degrees  K 

"c  \b  =  CHX  air  exit  temperature  m  degrees  K 

'e  al.a2.a3,bl,b2,b3  =  Coetficients  of  heat  capacities  tor 
"c  chamber  constituents  |02.N2,H20) 

%  The  simulation  ot  the  multivariable,  non  linear  equations  is  ill  be 
%  set  up  as  follows: 

% 

%  xdot=  F(x)*xdot  +  Gix.u)  -e  Hfx.wi; 

% 

%  Solving  foi  xdot  ; 

% 

%  \dot=  invil-Fix))*G(x.u)  ■<-  inv(I-F(x))*H(x,w) 


%  Setting  the  initial  conditions  of  the  .sutes  and  control  vanables 

x  =  134.464;  0.63420;2 1  1:0  0237911; 
u  =  j0.0025;0.3856;1.72e-3;1.5;t000|; 

% - 

%  Setting  the  other  variables 

Vc  =  157; 
al  =  0.9324; 
a2  =  0.90468; 
a3  =  1,831; 
bl  =  0.000184; 
b2  =  0.000204; 
b3  =  0.000744; 
qm  =  8586.33; 

R  =  8.314; 
xl  =  0.2344; 
vap  =  2287.0; 
hv  =  2213.0; 
mleak  =  0,0025; 


% 


%  Overall  number  of  integration  time  steps 

n  =  1200; 

p  =  6; 

r  =  p-1; 

%  T  is  the  sampling  duration  required  in  the  integration  of  the  equations 
%  of  motion  that  define  the  dynamics  of  the  processes  involved  in 
%  regulating  temperature,  humidity,  pressure  and  composition. 

T  =  0.004167; 
time  =  i  1  :n*p+  1  j; 
time  =  time»T; 
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%  Simulaiion  tor  n  integmtion  steps 

tor  k  =  in: 

:t  k-:  <  0 
t  =  0.0; 

'  delete  time.dat;* 
tprintt't 'time.dat'.  ’'re',  t) 
run  Itmn.inatUbIbrucej 

end; 


The  simulation  will  run  tor  p  integration  time  steps  prior  to  sampling 
%  and  evaluating  the  data  tor  parameter  estimation  and  control  input. 

for  i  =  p*k-r:p*lc 

02  =  0.3856; 
qh  =  2213; 

H20  =  0.67; 
x5  =  xi3,i)  -i-  273; 
x6  =  ui4.i)  f  281  33; 
x7  =  ui4,i)  +  8.33; 

xl  =  x(  1,1); 
x2  =  x(2.i); 

-i  =  x(3.i); 
x4  =  x(4.i); 
ul  =  ud.t); 
u2  =  u(2,i); 
u3  =  u(3,i); 
u4  =  u(4,i); 
u5  =  u(5,i); 

%  Calculate  the  partial  pressure  ot  H20  in  the  exit  stream  of  the 
%  CHX  and  the  accompanying  heat  of  vaporization 

Pswr  =  10'f28.5905 -8.2/2. 306*log(x6)  +  0.002484»x6-  3l42/x6); 

%  Calculate  the  heat  of  vaporization  of  the  moisture  in  the  exit 
%  stream  of  the  CHX 

hvapl  =  <13.11  -  log(Psur*0.9869))*R/i8.02'‘x6; 

if  x3  <  x7 
Pstar  =  0.0; 
end; 

%  Calculate  the  humidity  of  the  CHX  exit  stream 

we  =  18.02/(4*x2  +  28)*(100*Psur/(xl-l00*Psur)); 

%  Calculate  the  enthalpies  of  the  entrance  and  exit  gases  of  the  CHX 

CpTl  =  x2‘(al»x3  +  hl»x3'2)  -I-  (l-x2)*{a2*x3  +  b2*x3*2>  ...; 

+  x4*(a3*x3  +  b3»x3'2); 

CpTi  =  x2»(al»x7  +  bl*x7*2)  +  {l-x2)*{a2*x7  -t-  b2»x7*2)  ...; 

+  we»(a3»x7  +  b3*x7'‘2); 
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%  Calculate  the  change  in  sensible  heal  removed  across  the  CHX 


delh  =  CpTZ  -  CpTl  . 


It  delh  0; 

Jelh  =  0  0. 
emit 


'i  Calculate  the  change  in  absolute  humuiity  across  me  CHX 

delw  =  we  -x4t 

If  delw.  >  0: 
delw  =  0.0: 

end; 

‘b  Calculate  the  latent  heat  removed  across  the  CHx 

dell  =  delw^hvapl: 

Jelh2o  =  'u5<  ( 1  —  x4))*'delw: 

%  Calculate  the  total  heat  removed  across  the  heat  excahnger 

Jelch.x  =  lu.'h  1  +  x4))*{delh -rdell); 

if  delchx  >  0 
delchx  =  0.0; 

end: 

%  Calculate  the  function  Q 

Q  =  ix2*(al  -  al»x3ix5  ^  Z*bl*x3  -  hl»x3"2;x5i  ^  il  -  x2'|*  ...; 
lil  -  a2*x3/x5  i-  2*h2*.x3  -  b2*x3'2/x5)  +  x4*f-hv/x5 
a3  •  a3*x3/x5  +  2*b3"‘x3  -  b3»x3‘2/x5)); 

%  The  following  is  the  mathematical  represenutions  of  the  matrices 

%  F.  G.  and  H, 

F  =  (0.  0.  xl/x5,  0;  -x2/xl.  0,  x2/x5.  x2/(l  +•  x4);  -(x2*  .... 
(al*x3  +  bl»x3*2f  +  tl  -  x2)*(a2*x3  +  b2*x3*2)  -c  x4*  ...; 
ihv  -i-  a3*x3  +  h3*x3*2))/(xl*Q).- (al*x3  *•  bl*x3*2- 
a2*x3  -  h2»x3*21/Q,  0.  -((-l/fl  *  x4))«(x2*(al*x3  -i- 
bl*x3*21  -r  i;|  -  x2)f(a2«x3  +  b2*x3*2))  -r  a3*x3  .... 

b3*x3'2  -  hv  -  tx4/fl  +  x4))*fhv  +  a3*.x3  +  b3*x3*2)VO;  -t 
■l-(xl*il  x4-  l,(x4  +  ID),  0.  I.(x5*(l/x4  I  txd  +  In).  01; 

G  =  |R»x5'i(4*x2  *  28  +  x4*18)*Vc)*(u2  -e  u3  -  ul  -c  delh2o  1; 

‘  \  -  x4)*R*x5/((4*x2  +  28  +  1 8*x4)»Vc*x  1  )*(u2 -  x2*  .... 
ul),  (1  a-  x4)*R*x5/((4*x2  +  28  -t-  18*x4)»Vc*xl*Q)’'delch.x; .. 
il  -c  x4t*R*x5/f(4*x2  +  28  +  18.02»x4)*xl*x4*Vc*f l.'x4.... 

+  l.'tx4  +  ll))*u5/fl  +  x4i»delw|; 


H  =  (R*x5/((4'*x2  +  28  +  18*x4)*Vc)*(-02  +  H20  +  mleak);  ..  . 
R*x5*tl  ^  x4)/((4*x2  +  28  +  18*x4)*Vc»xl)*fxl*mleak-  02); 
R»(l  c.  x4)*x5/i(4*'x2  +  28  +  18*x4)«Vc*xl*0)*fqh  +  qm  + 
H20*vap);  R*x.5‘'‘f  1  4-  x4)/((4*x2  +  28  +  1 8»x4j*Vc*xl ‘.xd* . , , 
(l/x4  -  l.  (x4  +  1)))»H20|; 

%  The  following  sequence  will  calculate  the  change  m  x  over  the 
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^  integriiiion  time  pemni 


\doi  ==  inv!tjye(4)-Fj*'G  Hi. 

Kr.A-^  I  ■  ~  Xi  A  ■<  -  Kjot*T. 
a<  .1  1  »  —  u<  Ai. 


tfnd 


r’ 


‘^c  The  ing  sends  the  current  value  ot  x  into  the  t'ortran  parameter 

‘I  estimation  and  adaptive  control  routine,  wnnen  by  Dr  Cory  Finn,  that  has 
‘i  been  adapted  tor  this  tor  this  particular  control  problem 

'  del  y  dat;* 
y  =  XI  :.i^  !  i; 
tprintfCy  dat'.  're',  yi  1)); 
lpnnttry.dat'.  '%e'.  yilll; 
tpnntt'l’y.dat'.  'Te'.  yi3il; 
fpnntt't’y.daf.  '%e'.  yi4l); 

“i  The  following  sends  the  time  step  currently  in  clfect  to  the  Fortran  program 
!  del  lime  dat;* 

t  =  i»T; 

tprinttVhme.dat'.'^e'.i); 

%  The  output  of'  the  routine  is  the  new  control  inputs  that  minimize  the 
^  between  the  estimated  next  value  of'  x  and  the  setpoint. 

Iron  (finn.matlab|bnjce3 
load  inp.dat 
u(;,i+  i)  =  inp; 

ui  :.i-r  1)  =  u(  :.i); 

"r  ifk  >  25 

%  u(:,i+i)  =  [0  0025;0.3856;1.72E-03;0;i;00|; 

%  end; 


end; 


% - - - 

!  copy  lfinn.matlablbrl.datlfinn.matlablbr_y.dat 
'  copy  [finn.matlablbr2.dailfinn.matlablbr_u.dat 
!  copy  (finn.matlablbr3.dat  [fmn.matlablbr_e.dat 

!  del  (finn.matlablbrl.dat;* 

!  del  Ifinn  mailab|br2.dal;* 

!  del  Ifmn.matlablbr3.dat;* 


b.  Fortran  Code  Listings 
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1 .  Main  Program 


C - - - . - . . 

c . . . - 

c 

c  bhuce:-PLJuN't:  effective  lo,  :o  02 

c 

c . . - . 

c . . 

INCLUDE  in.NN  \nMO|COMMON  MIMO 

C  - . - . - . 

c 

C  INITIAUZATION  OF  v  ariables 

C 

c . - . 

KN  =  4 
KM  =  5 

INCLUDE  [FINN. PLANTIRAND  PLANT' 

OPEN(20.FIL£=(nNN.MATLABlTIME.DAr.STATUS=  OLD'.SHAREDI 
READCO.*)  TIME 
CLOSECO) 

IFiTIME.EQ.O.OlTHEN 
ITER  =  1 

OPEN(20.nLE-  |I=INN,MATLABlSTORE_SS2. DAT. STATUS-  UNKNOWN') 
READOO.*)  ITER.((Ya.;)J=«I.HDI).I-l.KN1. 

@  i(U(I.;).J=l.HDI),I  =  l.KM), 

®  ((THETAa..n.J=I.HDl),I-I.ICN). 

@  l(PHia,.r).F“I.HDl).I-l.KN). 

@  i((UMATRIX(U.IO,K  =  I.HDI),J=l.HDll.I  =  l.KN). 

@  .iDMATRIXa..n.J=I.HDl).I  =  I.KN). 

@  LAMDAlD.i-l.KNI 

CLOSECO) 

C  DO50J=I,HDl 

C  YiI.I)  =  J44W/I0.0 

C  YI2.D  =  0,<>4937*10.0 

C  YO.Ji  =  :i.i/io.o 

C  Y(4.I)  =  0.023279*500.0 

C  Ul  I. J)  =  0  0025*1000  0 

C  U(2.n  =  0.3856*10.0 

C  U(3.J)  =  0.00172*1000.0 

C  U(4,D  =  I.5*1.0 

C  U(5.F)  =  1000,0/1000.0 

C  50  CONTINUE 

C  DO  65  1  =  1.  KN 
C  DO50J  =  l.HDl 

C  D0  58  K=1,HD1 

C  UMATRIX(I.J.K)  =  0.0 

C  58  CON'HNUE 

C  DMATRIX(I.J)  =  1000000.0 

C  UMATRIXfl.J.J)  =  i.O 

C  PHKLJ)  =  0  0 

C  THETAd..!)  =0.1 
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C  60  CONTINUE 


C  65 

CONTINUE 

C' 

LAMDAUl  = 

1  0 

c 

L4MDA(:i  = 

1  a 

c 

LAMDAiii  = 

1  0 

c 

LVv1DA(4l  = 

;  0 

C  GOTO  1 00 
ELSE 

OPEN(;0.nLE=  |FINN  .VUTLAB|STORE.DAr.STATL'S«  UNKNOWN' 
READCO.*)  ITER.(( Y(I. J).J  =  1  ,HD1  i.l  =  1  KN>. 

■a  i(U(i.j).j=i.HDn.i  =  i.KMi. 

S  iiTHETA(l.DJ=I.HDn,I  =  l,KNV 

3  ':PHl(I,J).J=l.HDn.l  =  l,KN). 

S  ''(UMATRIX(I.J.Kl,K-l.HDn.J  =  l,HDli.i=  1  KN'i. 

@  'iDMATRIX(!.;)J-I.HDli.l=l.KN). 

9  'LAMDA(Ii.I=  I  .KNt 

CLDSEi:0l 

ENDIF 

C .  ... 

c 

C  SET  PROCESS  AND  REFERENCE-MODEL  PARAMETERS 
C 


NA(1)  =  1 
MB(1)  -  t 

NMOD(I)  *  NA(i;  • 

NAC)  =  5 
MB(2)  =  6 

NMODC)  -  NAI2)  -r  KM'MBCl 


NA(3)  =  5 
MB<3>  -  6 

NMOD(3)  =  NA(3)  -t-  KM*MB(3) 

NA(4>  =.  5 
MB(4)  -  6 

NM0D(4)  =  NA(4)  -i-  KM*MB(4) 

LAMDAMIN  =■  O.OOOOOI 
MEMORYLENGTH  =  1 .0 


sampletime  =  0.023 


C  TTNITin  =  l.2*NMOD<l)  •SAMPLETIME 

C  TINITt2)  -  1.2*NMOD(2)*SAMPLEnME 

C  TINITIS)  =  1 .2»NMOD(3)*SAMPLEnME 

C  TINITIA)  =  1  2*NM0D{4)*SAMPLETTME 

TPRED  =  6 

YHIn)  =  103.0/10,0 
YLO(l)  =  34.0/10.0 

YHJ(2)  =  0.64*10.0 
YLO<2)  =  0. 22*10.0 

YH1(3)  =  29.44/10.0 
YLO(3)  «.  15,53/10.0 
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YHI(4>  =  0.5;*500.0 

YIjO(4l  =  i)03975;*500,0 

UHlill  «  lUO'lOOO.O 
L'iX)in  =  0  0*1000  0 

USTEPll)  =  10*1000  0 

L'HICl  =  10  0*10.0 
UUOO  =  00*10.0 

USTEPCl  =  i  0*10.0 

UHllOl  =  10  0*1000. 0 
CLOiii  =  0.0*1000.0 

i;STEPi31  =  1.0*1000  0 

LHii4i  =  0*1,0 

ULX)l4i  =  0. 0*1.0 

USTEPI4I  =  0*1.0 

UHKS)  =  0500.0  1000.0 
UL0.51  =  0,0  1000.0 

USTEPfS)  =  OOO.O.'IOOO.O 


SETfll  =  34.464/10.0 
SETlOl  =  0  64937»10.Q 
SET(3i  =  Oil,  10.0 
SET(4t  =  0.003079*500.0 

[FiTIME.GE.SOlTHEN 
SET(l)  =  51.69/10.0 
SETIO)  =  0.437851*10,0 
SETO)  =  01.1/10,0 
SETf4)  =  0, 015337*500.0 
ENDIF 


C 

C  PLANT  OUTPUT.  PARANIETER  ESTftWTION  AND  C  ONTROL 

C 


100  CONTINUE 

OPEN!  00.  FI  LE  =  ■[  FINN .  MaTLABIY  .  DAT- .  STATUS  =  ■  UNKNOW .  SHARED! 

R£AD(20.*I  Yd. II.  YiO.ll.  Y(3,li,  Y<4,1| 

CUOSEIOO) 


Yd.l)  =  Yd, 11, 10. 0 
YlO.ll  =  Y(0. 11*10.0 
Y(3.1I  =  Y(3.U.'10.0 
Y(4.l)  =  Y(4. 1 1*300.0 

C  IF  (TIME.LE.TINITl  I ))  THEN 

C  Ud.l)  =  (0.0025-t-NRANDnTER+10)*0.0010)*lOOO  0 

C  U(O.l)  =  (0,3856-rNRAND(ITER+00)*0,0010)*10,0 

C  UrS.l)  =  (0.00170  +  NRAND(ITER+30>*0,0010)*1000  0 

C  U(4.1|  =  d.5  +  NRANDaTER+40l*0.010)*1.0 

C  U(5.1)  =  dO0O.0+NRANDfITER-*50)*I0.O)/lOOO.O 

C  ENDIF 

DO  1101  =  1. NAdl 
PHId.II  =  Yd.NAdl  +  O-II 
110  CONTINUE 

DO  113  I  =  1.NA(2) 

PHKO.D  -  Yf2.NA(2)  +  2-tl 
115  CONTINUE 
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DO  i:o  1  =  1  ,NA(3) 

PHi(3.I>  =  V(3.NAi3)  +  :-Ii 
'.30  CONTIN'UE 

DO  i:i  1  =  I.NAUi 
PHIi4J)  =  V'4.NA(4i*3-Ii 
i:i  continue 


DO  135  K=1.KN 
DO  i:;  i  =  i.MB(K) 

PHKK.l-NAlK))  =  Uil.MB(K)*3-II 
PH](K.I  +  NA(Ki  +  MB(IC))  =  U(3.MB<Kl-3-li 
PHl(K.l  +  NA(IO  +  3*MB(K))  =  U(3.MB<Ki*:i) 
PH](K.KNA<Kl  +  3*MB(K3)  =  U(4.MB(Kl-^3-li 
PHI(K.I4-NA(K)'^4*N{B(K))  =  U(5.MB(K)-:  II 
132  CONTINUE 
125  CONTINUE 

CALLVfULT(l) 

PREDERRORdl  =  Yil.li  -  CC 

CAU.  MULT  2) 

PRHDERROR(2)  =  Yi2.1i  -  CC 

call  MULTf3) 

PR£DERROR(3)  =  Y(3.1l  -  CC 

CALL  MULT(4) 

PREDERRORIA)  =  Y(4.U  -  CC 

CALL  UDU(  U 
CALL  UDU(2) 

CALL  UDU(3) 

CALL  UDU(4) 

C  130  IF  (TtME.LE,TtNIT(l)IC30TO  140 

135  DO  138  1  =  1,HD3 
DO  137J=1,HD3 
DIOPH(I.J)  =  0.0 

137  CONTINUE 

138  continue 

CALL  DIOPHANTINEI I .  U 
CALL  D10PHANTINE<1.2) 

CALL  DIOPHANTINE<1.31 
call  DIOPHANTINBT  .4) 

CALL  DIOPHANUNEd.S) 

CALL  D10PHANT1NE12.1) 

CALL  DIOPHANTINB2.2) 

CALL  D10PHANTINE(2.3) 

CALL  DIOKIANTlNEI2,4) 

CALL  D10PHANTINH2.5) 

call  DIOPHANTWEO.l) 
call  DiOPHANTlNEI3,2) 

CALL  DIOPHANnNEI3.3) 

CALL  DIOPHANTINE(3.4) 

CALL  DI0PHANT1NE(3.5) 

CALL  DIOPHANTINE(4.n 
CALL  DIOPHANTINE(4.2l 
call  DI0PHANTINE(4.3> 

CALL  DIOPHANTlNEI4,4) 

CALL  DIOPHANTINE(4.5) 
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CALL  TABLAU 


'Jd.n  =  CNEWMI 
=  L'NEWC) 

U(4.11  =  !;NW4i 
Ui5,n  =  i;neu-(3i 

!40  CONTINUE 

Call  display 

DO  150  1  =  1  KN 
DO  145  J  =  HDl,;.  l 
Yil./)  =■  Yd.Mi 
145  CONTINUE 
150  CONTINUE 

DO  160  1  =  1. KM 
DO  155  J  =  HD1.:.-1 
U(I.J)  =  UlI,J-l) 

155  CONTINUE 
160  CONTINUE 

ITER  =  ITER  *•  I 
■RITER.CE.OOOI  ITER  =  I 

190  OPEN(:0.RLE=  [FII-fN.NUTLAB|STCRE.DAr.STATUS=  unknown-) 
WRITE(:0.*)ITER.(  Y(l.J)J  =  I.HDl),I  =  l.KN). 

@  i(un.;).j=i.HDi).i=i.icM). 

@  i(THETAfI./).:=l.HDl).I  =  l.KN). 

a  i(PHia./).J=l.HDl).I=l.KN). 

a  (((UMATRIX(l.LK).K=l.HDn.J=l.HDl).l  =  l.K?n. 

a  !(DMATRIXn./l.J=l,HDl).I=I.KNI. 

a  (LAMDA(I).I=I.KNI 

CLOSEI20) 

OPEN(20.nLE='(FTNN.NUTLABHNP.DAr.STATUS=UNKNOWN-.SHAREDl 

wrjtei;o.»)  ud.U/iooo.o 

WRITEIOO.*)  U(2.l)/10.0 
WRITEI20.’;  UlJ.ll/1000.0 
WRITECO.*)  U(4.1)/1.0 
WRITECO.*)  Ul5,l)*l000.0 
CLOSE(20) 


009  END 


C - - - 

c— - - 


SUBROUTINE  DISPLAY 


INCLUDE  -inNN.MIMOICOMMON.MIMO- 


C  PRINTl.TIME.  Y(l.l)*10.0.  Y(2.1)/10.0.  Y(3.1)*10.0. 

C  a  Y(4.1)/500.0 

C  PRINT2.  Un.l)/1000.0.U(2.I)/IO.O.  Uf3.U/1000.0,U(4.I),  l  0. 
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c  a  U(3.i)*iooo,o 

C  PRINTS.  PREDERRORii  !  PREDERROR(:i.  PREDERRORii),  ?R£DERRORi4i 

!  FORMAT'  ,  4Fi:. 4.  Fi:  bi 
:  FORNLATT  5Fn.4i 
J  FORMAT'  ■,  4Fi:.4i 


OPEN(:O.P.LE=  IFINNMATLABlBRl.D.AT  .STAri.S=  ''SICNOWN' 
a  ACCESS  =  APPEND) 

OPEN(;i.nLE=  (FINN  MATLAB)BR:!  DAr.STATUS=  LNKNOVAN" 
a  -ACCESS  =  APPEND) 

OPEN(22.FILE='[nNNMATLABIBR3.DAT'.STATUS=  CNKNOW'N-, 
s  ACCESS  =  APPEND) 

WRITE(20.10)  TIME.  Vi  1.1)*  10  0.  Yi2.ii- 10  0.  V(3. 11*10.0. 

®  Yi4.1), 500.0 

WRITE(2l.ll)  TIME.  Ui  1.11  1000.0. U(2. 11.  lOO.UiJ.i)  1000  0. 

3  U(4.n.l-0.U(3. 11*1000.0 

WRITEl22.UmME.  PREDERRORl  1).  PREDERROR(2).  PREDERROR) Ji. 
®  PREDERRORl  4) 

10  FORMAT'  •.  4F12.4.FI2  6) 

11  FORMAT'  '.  4F10.4.FI4  4) 

13  FORMAT'  '.  5E14.4) 

CU3SE<20) 

CL0SB21) 

CIX)SE122) 

RETURN 

END 


C- 

c- 


INCLUDE  '(FINN.MIMOIMULT.MIMO' 
LNCLUDE  '(FINN.MIMOIUDU.MIMO' 

INCLUDE  '(FINN.MIMO|DIOPH_CPLD.MIMO' 
INCLUDE  (FINN.MIMOITABL_BRUCE2.MIMO' 
INCLUDE  [FINN.MIMOISIMPLEX.MIMO' 

C . . . . 

c . - . . . 


2.Subrouuoe9 


C 

C- 

C 

c 

c 

C- 

c- 


MULT.MIM02  EFFECTIVE  8/17/S9 


SUBROUTNE  MULTKKi 


C 
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uu  *  uuuouuu  ou 


c 


;  N  C  LU  D  E  '  1  n  N  N  NUMO I  CO  MMON  MIMO  ■ 


INTEGER  KX 


cc  =  0  c 


DO  350  1-1. NMOD(KK) 

CC  =  CC  -  PHKKK.Il'THET.^lKK.n 
350  CONTINUE 

RETL'RN 

END 


UDU.NIIMO:  EFFECTTV'E  8/17/89 


SUBROUTINE  UDU(KK) 


INCLUDE  rFlNN.MIMOICOMMON.MIMO' 


DOUBLE  PREaSION  F(HDl).  G(HDl).  BETA(HDl) 

DOUBLE  PREaSION  NU(HDI).  MU(HDI),  LGAIN/HDl).  ENTRY 

INTEGER  KK.  ROW.  COLUMN 


C 


DO  405  ROW=  1  .NMOD(KK) 

ENTRY  =  0.0 

DO  400  COLUMN-  1  .NMODfKK) 

ENTRY  =  ENTRY  -  UMATRD(rKK.COLUMN.ROWl*PHl(KK.COLUMN) 
400  CONTINUE 

FfROW)  -  ENTRY 
405  CONTINUE 

DO  410  ROW-  i  .NMODdOO 
G(ROW3  -  OMATRIX(XK.ROW)»F(ROW1 
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410  CONTINUE 


TRD  =  00 
SETAdi  =  1.0 

D0  415  J  =  l  .NMOD(KK) 

BETA(J  +  1)  =  BETA(J)  - 
TRD  =  TRD  +  DMATRIX(KK.J) 
415  CONTINUE 


MEMORYLENGTH  =  TRD 

LAMDAdOO  =  MEMORYLENGTH<MEMORYL£NGTH-^iPREDERROR(KKi**? 
a  BETAfNMOD(KKl+ll)) 

[RLAMDAdOO  GT.l.Ol  LAMDA<KK)  =  10 
IFiLAMDAlKKl.LT  LAMDANAN)  LAMDAdCK)  =  LAMDANUN 


BETA(l)  =  LAMDAlKK) 

D0  425  J=l.NMOD(KK) 

BETA(J-H)  =  BETA(J) 

DMATRIX(KK.T)  =  BETA(D*DMATRIXfKK.Jl,(BETA(J-*- IfLAMDAcKK.l 

NU(D  =  0<D 

MUlD  =  -RD.'BETAlD 

IF(J.EO  n  CXJTO  425 

DO  420  1  =  1, M 

UMATNEWd.n  =  UMATRlXlKK.l.n  -  NU(Il*MU(J) 

NU(I)  =  NU(I)  +  UMATRIX(XK.I.J)*NU(J> 

420  CONTINUE 
425  CONTINUE 

DO  435  Ul  .NMODlKKI 
LGAINII)  =  NU(I)/BETAlNMODlKKl-^I) 

DO430  J  =  l.NMOD(KIO 
UMATRDOKK.I.D  »  UNlATNEWd.D 
IF(l.EO  J)UMATRIX(KK,I.D  »  1,0 
IFO.LT.I)  UMATRD((1CK.I.D  •  0.0 
430  CONTINUE 
435  CONTWUE 

D0  440I  =  I.NMOD(KK) 

THETA(KK.i)  =  THET.AlKK.D  +  LGAlN(r)*PREDERROR(KK) 

440  CONTINUE 


RETURN 

END 


C . . . . . 

C 

C  DI0PH_CPLD.MIM02  effective  8/17, ’89 

C 

C . . . . 

C . . . 


SUBROUTINE  DIOPHANTINBKY.KU) 


C- 


98 


n  n 


INCLL  DE  I  FINN  MIMOirOMMON  M!MO 


WJL'BLE  PRECISiON  THiHDIi,  SiHDi.HDJ).  CONST^VNT.KDJi 
INTEGER  kV.  kU.  S!.  S,  DMAX.  DMIN 


INITTAUZE  THE  VARIABLES 


M  =  NAiKV! 
N  =  MBcKY/ 


DNUN  =  1 
OMA.X  =  NtBfKYi 


DO  510  1=  I.HDj 
DO  300J  =  I  HD3 
SiI.Jl  =  0.0 
500  CONTINUE 

CONSTANTin  =  0.0 
510  CONTTNUE 


tX)5:0I=l.M 
TH(I)  =  THETAiKY.l) 
5:0  CONTINL'E 


DO  5:i  I  =  M  +  l.Mi-N 
TH(I)  =  THETAiKY.N*(KU-I(-1i 
3:i  CONTINUE 

THlM*N*ll  -  0  0 


C . - . . 

c 

c  CALCULATE  THE  ALPHA  S  AND  BETA'S  CORRESPONDLNG  TO 
C  THE  past  OUTPUTS.  PAST  INPUTS.  AND  FUTURE  INPUTS 

f'  FOR  THE  PREDICTED  Y'S  FROM  Y(t+ II  TO  Y(t+r> 


DO  530  J  =  l.N-i-M 
Sil.Jl  =  TH(J) 

530  CONTINUE 

C  Sil.(M+N  +  TPRED+l))  =  THIM+N  +  U 
Sn.(M->-N*TPRED-M)(  =  0.0 


DO  580  1  =  2. M 

DO  550  J  =  1  .<M  +  N +TPRED) 

DO  S40K=1.(I-I) 

StI.J)  =  Sfl.JI  +  THIM-I  +  K+l)*SfK.T) 
540  CONTINUE 
550  CONTINUE 

DO  560J  =  I.M 
Sn.D  =  S(I.J)  +  TH(J-I  +  1) 

560  CONTINUE 

DO  570  J  =  il  +  M).{l4-.M  +  N-I) 
sn,j)  =  SILT)  +  TH(M+I) 

570  CONTINUE 
580  CONTINUE 
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DO  6:o:=<m-^iijTpred*i. 

D(3(>00J  =  l,'M-N-r?REDi 
DO  500  K»il  Mi.iM) 
i.I  Ji  =  Sil.il  *  TMiM-l  *K- !''S!K..'i 
rONTINUE 
CONTINUE 

0(J  510  J=.l  *Mi  i|  -  M*N  1, 

ST.Jp  =  S-I  Ji  *  THiJT»lp 
51  >3  CONTINUE 

5:0  CONTTNL'E 


DO  5:4  !  =  :.M 

do6::j=i..mi 

SiT,iM-N*-TPR£D-lii  =  SH.iM-^N  *  TPRED- I  h 
3  -  THai*S((I-Ji,iM'N-TPRED-l).i 

5;:  CONTINUE 

SlI.(M-N-TPR£D*l)‘  =  S<I.iM*N-TPRED-li.-TH'M-N-'.< 
5:4  CONTINUE 

DO  518  l  =  iM-<-  ImTPRED*  h 
D0  5:6J=I  M 

Sa.tM-i-N-TPRED»  I II  =  Si|.iM-N-TPR£D*  hi 
It  *  thi;)*S('I-;i,!M*n-tpred*  111 

5:6  continue 

Si|.i,M-^N»TPRED»  hi  =  S(I.<M-N  -  TPRED-  i  ii  ♦  Tri.  M  -  N  -  I . 
518  CONTINUE 

C . 

c 

C  CALCULATE  CONSTANTS  rOR  THE  CONSTRAINT  EONS 
C 

C . 


IFiKU  GE.:)GOTOS80 


DO  6501  =  1. iTPRED+l) 

DO6J0J=I.M 

CONSTANTdl  =  CONSTANTH)  ^  Sa.Ii'PHKKV./i 
6i0  CONTINUE 

DO  640  J » I M -p- 1  >.(  M  -  DNIAXI 

CONSTANTdl  =  CONSTANTiI)  *  Sil.Ji’PHKKY.Jl 
640  CONTINUE 

CONSTANTII)  =  CONSTANTd)  =  Sd.lM*N -TPRED- i  ii 
650  CONTINUE 

C)0  670I  =  I.TPRED 

DO  6601  =  1.  TP  RED 

DIOPH(TPR£D*<KY-n->-I.Ti  =  Sd*  l.J>^.M-DMAX) 

660  CONTINUE 

DIOPHITPRED'IKY  l)'=I.KM*TPRED+l)  =  CONSTANTd*  1) 
6'0  CONTINUE 

(iOTO  730 
680  CONTINUE 


DO  700t  =  l.(TPRED+li 

DO  OW  J  -  <  M  +  1 )  .1 M  +  DMAX) 

CONSTANTd)  =  CON'TANTd)  +  Sd.J)*PHl(KY.lKU-h*N  *li 
600  CONTINUE 
'00  CONTINUE 

DO  710  l-l.TPRED 
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DO  -SOJ  =  i  rPRED 

D!OPHiTPRED*'KY  1 1  - 1  J -.KU-i  .•TPREDi  =  vi  -  i 
roNTlNE’E 

DIUPH'TPRED'iKY  li  •  l,kM*TPRED-  I  -  ^ 

«  D10PH'TPRED"KY  ;  -  l.KM'TPRED  -  i  ■  -  '.'ONSTaN'!-; 

:'J  CONTINUE 


■JO  RETURN 
END 


C". 


c  . . . — . . . 

c 

c  tabl^brl'ce:  MiMo:  effective  lo o: 

r 

r . . . .  . 

r . - . . . . 


SUBROUTINE  TABLaU 


INCLUDE '(FINN-NUMOICONEMON, SUMO' 


INTEGER  KY.  KU,  M.  N,  T,  TAG 


T  =  TPRED 

TAG  =  1 

N  =  KN 

M  =  tCM 

SOO  !ER  =  •) 

DO  SIO  1  =  1. HD4 
DO  805  J=I.HD5 
Ail.Ji  =  0.0 
805  CONTINUE 
BlO  =  0 
BORJG(I)  =  0 
Dill  =  0  0 
W(l)  =  0  0 
810  CONTINUE 

DO  815  l  =  l  .HD5 
CiF)  =  0.0 
815  CONTINUE 


DO  830  KU  =  i  .M 

.Ai(KU-l)*T^  l.<KU-l)*0*TJ-i)=  1  0 
AcKU-lUT-t-  l  .iKU-h*I*T*T*n-  1.0 
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C  KL  i  *:»T-!  =  -I  STEP' Kf  I  -  L'Kl. 
c.  KI-;  TSTEP-KUi  : 

T 

\  '.  KU  ii-:*?-'.*  ,  0 

s-kU-ii-T-:  ;  !■>  =  -t  o 

s-  KU  i  I'T- I  KL'  ;  i-:"T  -  r- (.  =  jO 

A  AU  i)’:<T*T-j.=  i  ^ 

=  i  STEP-Kf' 

r  KT  s  steP'KUi 

5:0  c'OSTiN’UE 
x>  iZ!  ;  =  I  T 

A :’Ki. . :i*T - ;  :'M‘T- :;*s*T-  kc  ;  •;•?-!'=  1  .• 

A  Ki -I:*!-;  ;’M’T*;*N*T-.Ki.' : 
CM*M*T-;*N*T*-kU-li":*T-J,=.  UU)'ICU. 

C  :*M*T-;*NT*  KU-I  THJ.KUi 

s:5  CONriNfE 
S30  CONTINUE 

00  845  KY  =  UN¬ 
DO  S40J=  IT 
DO  555  [  =  ;  M‘T 

A.I,;»M*T-.K5-  =  T.AG*DIOPHt'KY  ^ 

A  I,:*M*T-'K;Y-b*:'T-T-J)=  TAG-DIOPHiiKY  h'T-J  I, 
A.1.4*M*T*:*N*T-'KY-n*:*T*Ji=  OiOPHfiKY.ii’T-i  ii 
A  I,4*M*T-:*N*T-  KY1'*:*T*T*J)=  DIOPHi'KV- b'T -  J  1, 
S35  CONTINUE 

C.:*M*T-'KY-ii':*T-Ji  =  TAGNYIjOiKY)  ■ 
a  OiOPHi(KY-li*T*J.M"T*  I'i 

C:*M*T*-KY-;i‘;*T-T-J..  TAC-f-YHIiKYi  - 
S  D10PHi(KYIi*T-J.M*T-i,i 

A(M*T-'KY-l)«:*T-J.4*M*T*:*N»TA'KY-r>*:*T-fi=  10 
AiM*T-(KY-U‘:*T-T-'J  4*M*T-2*N*T»<KY1i*;*T'A“  10 
A(M»T  +  (Ky-l)*:*T-J.4*M*T*:»N*T^(KY-I>*:*T*T-Jl'»  -10 
AiM*T-(KY-Ii*:«T-'-T-'-J.4»M*T*:*N*T*.KY-ll*:*T*T-;)»  l  0 
04*M»T-'-:*N»T*(KY-1)*2»T+«-  SETHCY)  - 
a  DIOPH((KY-I)*t-'-J.M»T-*lt 

C(4«M*T*:*N*T-iKY-l»*0*T^T-n»  -SET(KY)  - 
8  DIOPH(OCy-l)*T-'-J.M*T-k|. 

AiM*T-KKy-|i*:’T-'-J,4*M*T-'-4*N*T-MKy-IT:*T-Ji=.  1  0 
aiM*t-(ky-1}*:*t-t-i.4*m*t-4*n*t*iky-ii':«t-t-;.=  i  o 
Ci4*M*T+4*N*T*<Ky-ll*2*T-'-T(»  0.0 
C(4•M•^-4•N•T-'Ky-l)•2•T■^T-^Il=  0  0 
840  CONTLNUE 
845  CONTINUE 

DO  850  1=  I  T 
D<I)  =  0  000001 

850  CONTINUE 

DO  851  I  =  iT~H,:'T 
Dd)  =  0.000001 

851  CONTINUE 


D0  85:1  =  i2*T-'-I).5*T 
txn  -  0  000001 

852  CONTINUE 

D0  853I«'3*T*li  4*T 
Ddl  •  u  OOOOOl 

853  CONTINTJE 

DO  854  l=i4*T*  1  |.5*T 
D<!)  »  0  '30000! 

854  CONTINUE 

DO  *55  I-1.2»N*T 
D(UM*T)  -  1,0 
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S55  CONTINUE 


(■  ADD  EXTRA  artificial  VARIABLES  FOR  PHASE  i  OF  THE 
C  SIMPLEX  ROUTINE.  .SONIMIZE  THE  SUM  OF  THESE  NEU 
C  variables  to  obtain  a  basic  feasible  SOLUTION  IF 
C  THE  OBJECTIVE  FUNCTION  IS  CREATER  TH.AN  ZERO  THEN 
C  NO  FEASIBLE  SOLUTION  EXISTS 
C 

c . - . . . 


DO  8601  =  1. iM»T*:*N*Tl 
Aa.A*M*T*6»N*T-U=  VO 
Ci-»'M*T-6»N»T*li  =  I  0E03 
Bill  =  4*M*T»6*N*T-I 
560  CONTINUE 


MCON  =  M"T-:*N’T 
NVAR  =  5*M*T-8*N*T 


call  simplex 

DO  875  1  =  1, MCON 
DO  870  1  =  1. MCON 

IF  (B(J).EQ  BORJG(Il)  THEN 
IFlDlJl.CT.l  0E-08)THEN 
lER  =  100 
PRINT* 

PRINT865.  lER 

86S  FORMATC  error  T6.':  No  fcaiiblc  »oluUons  i 

ENDIF 
ENDIF 

870  CONTINUE 
875  CONTINUE 

IF  (1ER.NE.0)  THEN 
IF  (Tag.eq.o)  then 
D0  87JKU«I.M 
UNEWIKUI  =  U<KU.;i 
878  CONTINUE 
GOTO  890 
ENDIF 
TAG  =  0 
goto  800 

ENDIF 

DO  885  KU  =  1.M 
DO  880  1=1. T 

UPRED(KU.I)  =  W(iKU-l)*T  +  ii 
880  CONTINUE 

UNEWfKU)  -  W((KU-I>*T+I) 

885  CONTINUE 

PRINT* 

PRINT*.  OBJ.  (UNEWCKUI.  KU=1.MI 

390  CONTINUE 

D0  894I  =  I,HD3 
D0  892J«I.HD3 
DIOPH(t.J)  =  0.0 
393  CONTINUE 
894  CONTINUE 


RETURN 
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o  o 


END 


C . . . 

C . . . 

c 

C  SIMPLEX.NUMO;  effective  8.  16  89 
C 

C . . . 

c - - - - 


SUBROUTINE  SIMPLEX 


LNCLUDE  ■[  RNN .  MIMO|COMMON .  NflMO 

DOUBLE  PREaSION  KMIN.  RMAX.  .^VHAT.  EPSl.  EPS? 
INTEGER  K.  R.  NMULTP.  NMULTD 


EPSl  -  l.OE-08 
EPS?  =  l.OE-W 

DO9001«l.MCON 
BORlOm  >  B(I) 

900  CONTINUE 

D0  910I  =  1.NVAR 
A(MC0N  +  1.I)  »  0,0 
DO905  J»l.MCON 

A(MCON+1.0  =  A(MCX)N+I.I)  +  C(80R10(T»)*A(J.I) 
905  CONTINUE 

aimcon + 1  ,n  =  A(McoN + 1 .1)  -  cm 

01 0  CONTINUE 

D(MCON  +  l)  =  0.0 
D0  91S  J-l.MCON 

D(MCON  +  I)=  D<MCON  +  l)+C(BORlG<D)*Dm 
915  CONTINUE 

NMULTD  «  0 
NMULTP  «  0 


920  KMIN  =  -EPS2 
K  ="  0 

D0  925I=I,NVAR 
IF  IA(MCDN+ 1  .n.LT.KMIN)  THEN 
KMIN  »  AlMCON  +  l.n 
K  =  I 
ENDIF 

925  CONTINUE 


IF  nC.EO  O)  GOTO  960 
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RMAX  =  1  0E30 
R  =■  0 

DO  “liS  1=  I  .MCON 

!F  ((Ad.Ki.GT.EPSlI  AND  iiDili.AiI.K>i  LI  RNtAMi 
li  THEN 

RMAX  =  D<n  Ail.K) 

a  -  1 

ELSEIF  KAiI.Ki  OT  EPSn.AND.(iD(Ii  Aii.K.ii  EQ.RMAX. 
S  THEN 

DO  030  J=  I.N'VAR 

IF  i(A(I.Ji.AiI.Kii,GT.iAiR.Ji.  AcR,K»)GOTO035 

!F  UAd.JV  Ail.K)).EO  i  A!R,H  AiR.KiiiGOTO  030 
RMAX  =  Ddi  AiI.Kj 
R  =  I 
GOTO  035 
030  CONTINUE 
ENDIF 

035  CONTINUE 

IF(R,EQ.O)  THEN 
lER  =  101 
PRINT* 

PRINT940.  lER 

OAO  FORMATC  ‘.'ERROR  ‘16.’:  The  solution  is  uoboundcii'i 
GOTO  960 
ENDIF 

DO  9501  =  I.MCON*! 

IF  (I.EQ.R)  GOTO  950 
AHAT  =  A(I.K>/A(R.K) 

IF(DASS(AHAT).LT.EPS1|AHAT  =  0.0 
Dfl)  =  D(I)-D(R)*AHAT 
IF(DABS(D<I)).LT.EPSI)D(I)  =  0.0 
D0  94S  J=I.NVAR 
AH.D  =  A(I.J)-(A(R..n*AHAT» 
IF(DABS(Aa.D).LT.EPSl)A(I.fl  «  0.0 
945  CONTINUE 
950  CONTINUE 


AHAT  =  A(R.I0 
D0  955J=1.NVAR 
A(R.D  =  A(R..n/AHAT 
IF(DABS(A(R.D).LT.EPS1IA(R.D  =  0.0 
955  CONTINUE 

DlRl  =  D(R)/AHAT 
[F(DABS(D(R)|.LT  EPSl)D(R)  =  0  0 

SIR)  =  K 

GOTO  920 

960  CONTINUE 


OBJ  =  D(MCON  +  n 

DO  9701  =  I.MCON 
WO)  =  0.0 
DO  965  J- I.MCON 
W(I)  =  Wd)  +  CIBOl'Ad.BORIGdU 
965  CONTTNUE 
070  CONTTNUE 


C 


DO  9S0  I -I.MCON 


C  IF  (DABSiDdli.GT  EPSIIGOTO  080 

C  DOo’g  K=  I  NVAR 

C  DQo-'s  j^i  mcoN 

('  ;F  K  EQ.BiJii  GOTO  o’g 

G  IF  iK  EO.BORlGlJu  GOTO  o-g 

C  0-6  CONTINUE 

r  !F  .D.ABSiAil.Kii  LT  EPShGOTO  0-8 

C  AHAT  =  AiMCON  +  I,K6Ai1,Ki 

C  DO  0’4J  =  I  NVAR 

C  IF  KAfMCON^  1.JHAHAT*Ai!.J)()  LT  EPS!) 

C  GOTO  9^8 

C  974  CONTINUE 

C  NMULTD  =  NNOILTD  -  1 

C  PRLNT* 

C  PRINT97:,  NNa.'LTD.  I.  Bin.  K 

C  972  FORNtATC  '.'Multiple  DuAi  Solutions  .4141 

C  978  CONTINUE 
C  980  CON'nNUE 
C 

C  DO990I=1.NVaR 

C  D0  985J  =  i.MC0N 

C  IF  II.EQ.BO)  GOTO  900 

C  985  CONTINUE 

C  IF  (DABSfA(MCON-l.Il).GT.EPSnGOTO  990 

C  NMULTP  =  NMULTP  +  I 

C  PRINT* 

C  PRINT98g.  NMULTP.  I 

C  988  FORMAT!'  '.'Multiple  Primal  Solutions'. 2I4> 

C  990  CONTUNUE 


RETURN 

END 


C - 
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APPENDIX  B.  PARAMETRIC  STUDIES 
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Coolant  Flowrate  through  Heat  Exchanger 
as  a  Function  of  Chamber  Pressure 


jH/33|  ib3H  M^noiqj  aiBJMOU  lUBpoo 
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Chamber  Pressure,  psia 


Coolant  Temperature  into  Heat  Exchanger 
as  a  Function  of  Chamber  Pressure 
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Chamber  Pressure,  psia 
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3.  ‘i3qureq3  ojui  axnieidduisx 
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Chamber  Pressure,  psia 


Absolute  Humidity 
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Chamber 


Absolute  Humidity  of  Air  out  of  Heat  Exchanger 


ire  /Cip  2y{/i9ve/t<  ‘isSueq^xg  resH  JO  jno  iiy  Jo  XjipiiunH  sqv 
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(.'harnber  Pressure,  psia 


Chamber  Absolute  Humidity 
as  a  Function  of  Chamber  Pressure 


ire  itip  iyvmtti  Sjj  ‘XiipnunH  amiosqy  i3qureq3 
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Chamber  Pressure,  psia 


Heat  Exchanger  Area 
as  a  Function  of  Chamber  Pressure 


jUi  ‘Baiv  wSucMOxg 
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Chamber  Pressure,  psia 
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